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ABSTRACT 

The  effects  of  moderate  nonlinearity  on  the  propagation 
of  sound  are  appreciable,  and  become  dominant  at  very  high 
amplitudes.  These  effects  and  the  phenomena  of  linear 
acoustics  are  described  by  the  second-order-nonlinear  wave 
equation,  which  is  derived  in  this  thesis  and  solved  by 
numerical  means.  The  validity  of  the  solution  is  demon¬ 
strated  by  its  agreement  with  various  approximations  in 
their  domains  of  applicability,  and  by  its  reproduction  of 
results  derived  from  experiments.  Using  the  numerical 
solution  in  simulation  of  the  operation  of  acoustic  trans¬ 
ducers  at  finite  amplitudes,  conclusions  are  presented 
concerning  the  amount  of  energy  that  may  be  transmitted  to 
the  far  field  by  various  types  of  signals. 
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according  to  the  following  convention: 

*  :  asterisk,  used  to  denote  a  reference  quantity; 

o  :  denotes  a  nominal  quantity,  or  one  referred  to 
a  mean-carrier  frequency; 

1  :  denotes  the  higher  carrier  of  a  bifrequency  pair; 

2  :  denotes  the  lower  carrier  of  a  bifrequency  pair; 

-  :  denotes  the  difference  frequency  of  a  bifrequency 

pair; 

n  :  denotes  the  n-th  component,  e.  g.,  of  a  spectrum. 

For  example,  the  variable  u  is  a  particle  velocity, uA 

a  reference  velocity,  u  a  nominal  or  mean-carrier  particle 
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velocity,  and  so  on.  In  certain  instances  a  subscripted 
symbol  has  a  different  meaning;  such  symbols  are  given 
explicitly  in  the  list  of  symbols  which  follows. 

The  superscript  prime  (')  is  used  to  denote  an  excess 
quantity  above  the  nominal  or  equilibrium  value,  unless 
otherwise  defined.  For  example,  the  excess  acoustic 
pressure  is  denoted  p'=p-p  . 
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Chapter  1 

INTRODUCTION 

Acoustical  signals  may  easily  exceed  those  levels  for 
which  the  linear  lossy  wave  equation  provides  an  adequate 
model.  Even  if  the  amplitude  is  small,  distortion  products 
will  accumulate  unless  they  are  removed  by  absorption  in 
the  medium.  The  wave  equation  may  be  taken  to  second  order 
to  include  the  effects  of  moderate  nonlinearity. 

Numerous  solutions  exist  for  particular  forms  of  the 
second-order-nonlinear  wave  equation,  for  example,  Burger's 
equation  (1948).  Analytic  solutions  have  been  obtained  only 
for  restricted  cases,  such  as  for  particular  regions  of  the 
field  or  for  plane-wave  propagation. 

The  purpose  of  this  thesis  is  to  present  an  original 
numerical  solution  of  the  second-order-nonlinear  acoustic 
wave  equation,  applicable  to  plane-wave  propagation  and  to 
propagation  from  an  axisymmetric  source  of  finite  extent 
having  an  arbitrary  amplitude  and  phase  profile.  As  the  solu¬ 
tion  is  a  numerical  procedure,  no  approximations  are  made  in 
the  second-order-nonlinear  acoustic  wave  equation  as  it  is 
solved.  This  represents  an  advance  on  the  prior  art,  in 
which  solutions  were  obtained  for  simpler  forms  of  the  equa¬ 
tion  solved  in  this  thesis.  The  numerical  solution  presented 
in  this  thesis  is  valid  within  an  unbounded  medium  which  is 
either  lossless  or  has  any  desired  attenuation  character- 


istic  as  a  function  of  frequency,  and  is  either  non-disper- 
sive,  or  has  any  desired  phase  velocity  as  a  function  of 
frequency.  The  useful  conditions  of  a  t he rmoviscous  or  a 
monorelaxing  medium  are  special  cases  which  may  be  handled 
with  ease. 

The  first  chapter  of  this  thesis  introduces  the  funda¬ 
mental  concepts  of  nonlinear  systems  and  presents  a  review 
of  the  prior  art  in  solutions.  A  brief  derivation  of  the 
second-order-nonlinear  wave  equation  is  presented  in 
Chapter  2  and  a  numerical  solution  of  this  equation  in  Chap¬ 
ter  3.  The  results  of  the  solution  under  various  conditions 
are  discussed  in  Chapter  4,  and  some  of  its  implications  are 
presented  in  Chapter  5. 

Nonlinear  Systems  and  Equations 

A  nonlinear  system  is  one  for  which  the  principle  of 
linear  superposition  fails.  That  is,  if  f (x)  is  the  output 
of  the  system  in  response  to  the  input  x,  then 

f(ax+by)  f  af(x)+bf(y)  (1.1) 

The  functional  relationship  between  two  variables — 
for  example,  the  pressure  on  and  condensation  of  an  ele¬ 
ment  of  fluid — may  be  expressed  as  a  power  series.  Ifp'  is 
the  excess  acoustic  pressure  and  s  the  condensation,  then 
the  series  in  Beyer's  notation  (i960)  truncated  at  the 
second  term  may  be  written 
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p*  =  As+(B/2)s2,  (1.2) 

which  is  called  second-order  as  the  first  missing  term 
is  proportional  to  the  third  power  of  s ,  a  number  much 
less  than  unity.  The  dependence  of  pressure  on  the 
specific  entropy  has  been  omitted  from  Equation  (1.2),  as 
the  latter  is  intended  to  apply  for  values  of  the  Mach 
number  e  less  than  1/10.  The  Mach  number  is  the  ratio  of 
the  peak  particle  velocity  to  the  small-signal  speed  of 
sound.  Whitham  (1974)  has  shown  that  the  change  in 
entropy  across  a  shock  front  is  of  the  order  of  the  Mach 
number  cubed,  and  is  therefore  negligible  in  a  second- 
order  analysis. 

The  subject  of  this  thesis  is  a  numerical  solution  of 
the  second-order-nonlinear  wave  equation  in  several  systems 
of  coordinates.  A  numerical  solution  is  a  procedure  that, 
given  the  values  of  the  dependent  variables  at  one  set  of 
values  of  the  independent  variables,  specifies  how  the 
former  may  be  computed  at  another  set  of  the  latter. 
Numerical  solutions  generally  involve  a  degree  of  error 
due  to  approximation,  which  may  be  reduced  as  far  as  is 
desired  by  appropriate  numerical  techniques. 

Previous  Solutions  to  Nonlinear 
Wave  Equations 

This  section  reviews  previous  solutions  to  nonlinear 
wave  equations.  Beginning  with  Euler's  equations,  Earnshaw 


(1860)  obtained  a  result  valid  for  inviscid  progressive 
plane-wave  propagation  and  determined  that  the  velocity  of 
propagation  is  a  function  of  amplitude: 

c=c+$u,  (1.3) 

o 

where  g  is  the  parameter  of  nonlinearity  of  the  fluid  after 
Beyer's  notation,  u  is  the  particle  velocity,  and  cq  is  the 
small  signal  speed  of  sound.  Riemann  (1860)  independently 
obtained  a  solution  which  includes  plane  waves  traveling  in 
two  opposite  directions. 

If  a  sinusoidal  signal  of  small  amplitude  (compared 
With  the  equilibrium  pressure)  is  followed  as  it  propagates 
through  an  inviscid  fluid,  the  compressions,  whose  particle 
velocities  are  positive,  will  travel  faster  than  the  rare¬ 
factions,  whose  particle  velocities  are  negative.  Thus 
adjacent  compressions  and  rarefactions  will  approach  one 
another,  and  a  discontinuity  in  the  pressure  will  be  formed. 
If  a  numerical  procedure  based  directly  on  Equation  (1.3)  is 
used,  at  the  point  at  which  a  discontinuity  is  formed  it  will 
fail,  as  the  compressions  and  rarefactions  continue  to  travel 
past  one  another,  and  the  predicted  pressure  will  have 
several  distinct  values  at  each  point  in  the  vicinity  of  the 
discontinuity.  This  method  is  capable  of  giving  good 
results  out  to  the  critical  range,  as  reported  by  Pestorius 
and  Blackstock.  (  197  3).  Beyond  the  critical  range  the 
authors  average  across  the  shock  front,  thus  suppressing  the 
multivalued  waveform. 


The  Fubini  (1935)  and  inviscid  Fay  (1931)  solutions 
for  initially  monofrequency  waves  give  analytic  solutions 
for  the  near  and  the  far  fields  of  the  lossless  Burger's 
equation.  In  the  near  field,  i.e.  for  0<1,  Fubini's 
expression  for  the  n-th  harmonic  of  the  pressure  is 

Pn(0)  =  if  Jn<na>’  0.4) 

and  in  the  far  field,  for  a>3,  Fay's  result  is 


P  (a)  = _ 

n  it  (1+0) 


(1.5) 


where  o=g£kz  is  the  distance  relative  to  the  plane-wave 
shock  formation  distance 


rc=  1/Bek. 


(1.6) 


The  Fubini-Fay  solution  for  the  near  field  and  the 
asymptotic  far-field  spectrum  of  a  decaying  shock  wave  were 
combined  by  Blackstock  (1966).  His  formulation  uses  a  poly¬ 
nomial  in  0  to  match  the  values  and  derivatives  of  the  near- 
and  far-field  solutions.  This  formula  gives  excellent 
agreement  with  comparatively  lossless  plane-wave  experimental 
data  throughout  the  field. 

The  limits  of  classical  theory  lie  in  its  neglect  of 
dissipation.  Fox  and  Wallace  (1954)  obtained  a  pertur¬ 
bation  solution  for  plane  waves  in  a  lossy  medium.  This 
solution  perturbs  the  rates  of  harmonic  generation  and  decay 
of  the  lossless  plane-wave  solution.  Their  paper  discusses 


the  usefulness  of  scaled  coordinates  such  as  will  be  used  in 


this  thesis  (though  their  notation  is  different),  and  also 
considers  the  impossibility  of  generating  acoustic  signals  of 
arbitrarily  large  amplitude  due  to  saturation  of  the  medium. 
Their  paper  also  reports  on  comparison  between  numerical 
simulation  and  experiments  in  air,  water,  and  carbon  tetra¬ 
chloride.  The  agreement  between  simulation  and  experimental 
results  is  excellent. 

Cook  (1962)  described  an  iterative  numerical  procedure 
for  the  calculation  of  the  distortion  of  plane  finite- 
amplitude  waves  in  a  lossy  medium.  This  procedure  is 
similar  to  that  used  in  this  thesis  for  nond ispers ive  plane- 
wave  propagation,  but  is  somewhat  different  in  its  details 
and  motivation.  The  model  is  based  on  two  assumptions:  that 
the  distortion  mechanism  may  be  described  by  a  change  in 
phase  velocity  which  is  directly  proportional  to  the 
particle  velocity,  and  that  the  absorption  of  each  frequency 
component  is  proportional  to  its  amplitude  times  the  square 
of  its  frequency.  The  procedure  is  to  distort  the  waveform 

as  it  propagates  over  a  small  interval,  and  then  to  correct 

; 

for  absorption. 

Another  nonlinear  wave  equation  of  great  interest  is 
the  Kor t eweg-DeV r ie s  equation;  as  given  by  Lamb  (1965),  it 
may  be  written  in  the  following  scaled  form: 

.82  _  p  9P  =  „  a3p 

do  3T  3T 3  (1.7) 

This  equation  describes  propagation  in  lossless  and 
dispersive  nonlinear  media.  As  Lamb's  analysis  (  19  6  5)  of 


this  equation  indicates,  the  signal 

P(t)  =  asech2 (t/D) ,  (1.8) 

where 

D  =  (12K/a)1/2  (1.9) 

is  a  steady-state  solution  of  the  K-dV  equation.  This  signal 
is  termed  a  solitary  wave  solution,  or  soliton.  It  will  pro 
pagate  through  the  medium  at  a  velocity  dependent  on  its 
amplitude  a.  As  Zabusky  (1967)  has  shown  by  means  of  a 
numerical  analysis  of  Equation  (1.7),  several  solitons  may 
interact  without  losing  their  identities.  Any  signal  pro¬ 
pagating  in  a  nonlinear  dispersive  medium  will  be  resolved 
into  one  or  more  solitons.  Under  certain  conditions  these 
may  coalesce  at  a  later  time  to  re-form  the  initial  signal. 

Rosen  (1966)  discussed  the  computational  solution  of 
nonlinear  parabolic  differential  equations  by  linear 
programming.  This  method  involves  choosing  a  set  of 
functions  of  which  a  linear  combination  approximates  the 
desired  solution.  Two  sets  of  conditions  need  to  be 
satisfied:  the  boundary  conditions  and  the  partial 
differential  equation.  One  may  choose  the  functions  so 
that  each  satisfies  the  boundary  conditions  and  a  linear 
combination  satisfies  the  differential  equation,  or  choose 
them  so  that  each  satisfies  the  differential  equation  and 
a  linear  combination  satisfies  the  boundary  conditions. 

In  either  case,  it  is  a  linear  programming  problem  to 
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determine  the  coefficients  of  the  linear  combination  which 
minimize  the  error  in  some  sense.  Of  particular  interest 
to  the  subject  of  this  thesis  is  that  the  method  is 
directly  applicable  to  the  plane-wave  form  of  Burger's 
equation.  Rosen  gives  several  examples  of  the  error 
bound  which  may  be  expected,  showing  that  with  as  few  as 
five  well-chosen  functions  an  error  of  less  than  1/10 
percent  may  be  obtained. 

The  work  discussed  above  dealt  with  the  propagation 
of  finite-amplitude  plane  waves  only.  It  is  possible  to 
use  a  stretched  coordinate  system  for  Burger's  equation 
which  will  accommodate  plane,  cylindrical,  or  spherical 
spreading.  This  approach  is  used  by  Fenlon  (1971)  in 
a  method  for  computing  the  interaction  between  spectral 
components  in  progressive  finite-amplitude  waves.  In  this 
method,  the  spectral  representation  is  truncated  at  a  finite 
number  of  terms,  and  a  coupled  set  of  nonlinear  differential 
equations  is  obtained  for  the  component  amplitudes.  Given 
an  initial  spectrum  at  a=o,  the  spectrum  at  o=Ao  may  be 
obtained;  then,  at  0=2Aa,  and  so  on  until  the  desired  range 
is  reached. 

Cary  used  the  transformation  of  Naugol'nykh  (1963)  to 
obtain  an  equation  from  which  an  expression  for  the  "extra 
decibel  loss"  due  to  finite  amplitude  absorption  was 
obtained  (1967).  This  paper  reaches  the  conclusion  that 
finite-amplitude  losses  in  spherical  waves  are  less  than  for 
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plane  waves  of  the  same  source  level,  but  still  are  not 
negligible. 

A  second  paper  by  Cary  (1968)  that  is  of  interest 
in  the  history  of  numerical  solutions  to  nonlinear  wave 
equations  presents  a  numerical  solution  based  on  Burger's 
equation  in  a  stretched  coordinate  system,  suitable  for 
plane,  cylindrical,  and  spherical  spreading.  Following 
the  method  of  Banta  (1965),  the  normalized  particle 
velocity  is  represented  as  a  Taylor's  series,  and  enough 
terms  are  kept  to  ensure  an  adequate  model  of  the  distortion 
process.  This  method  is  valid  until  the  effects  of  absorp¬ 
tion  have  become  dominant  over  nonlinear  effects.  Reason¬ 
ably  good  agreement  between  experimental  data  published  by 
Romanenko  (1959)  and  the  numerical  predictions  is  shown. 

Cary  notes  that  this  method  should  be  useful  in  the  design 
of  nonlinear  acoustic  systems. 

The  pressure  field  radiated  by  an  acoustic  source  may 
be  conceptually  divided  into  three  zones.  Zone  I  is  the 
region  close  to  the  projector  in  which  energy  is  transferred 
into  harmonic  or  modulation  frequencies  by  nonlinearity 
faster  than  it  is  removed  by  viscous  losses  or  reduced  by 
spreading  losses.  Conversely,  zone  III  is  the  region  far 
from  the  projector  in  which  either  viscous  losses  or 
spread ’ ng  losses  are  dominant  over  the  effects  of 
nonlinearity.  Between  zones  I  and  III,  a  shock  wave  may 
form  and,  if  so,  the  rates  of  harmonic  generation  and  loss 
will  be  comparable  within  some  region,  giving  rise  to  a 
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quasi-stable  waveform.  This  region,  if  it  exists,  is 
called  zone  II.  Cary  (1973)  has  published  an  exact  zone  II 
solution  of  Burger's  equation  for  a  parametric  source. 

Sadchev  and  Seebass  (1973)  published  a  study  of  the 
decay  of  spherical  and  cylindrical  shock  waves.  This 
article  presents  a  numerical  finite-difference  solution  to 
a  general  form  of  Burger's  equation  including  the  effects 
of  viscous  absorption  in  either  plane,  cylindrical,  or 
spherical  geometry.  The  plane  wave  form  may  be  transformed 
into  a  linear  equation  by  the  methods  of  Hopf  (1950)  and 
Cole  (1951).  The  latter  equation  has  been  solved  by 
Lighthill  (1956),  and  Sadchev  and  Seebass  use  Lightill's 
analytic  results  as  a  test  for  their  numerical  method. 

The  numerical  method  of  Sadchev  and  Seebass  employs  a 
mesh  of  sample  points  in  two  coordinates,  which  may  be 
considered  as  along  and  normal  to  the  direction  of 
propagation.  Cylindrical  and  spherical  spreading  each  have 
t  h  .■  effect  of  transferring  energy  outward  from  the  beam 
axis.  For  this  reason  it  is  necessary  to  increase  the 
number  of  mesh  points  in  the  normal  direction  and  the  size 
of  the  normal  mesh-point  spacing  from  time  to  time,  so  as  to 
represent  the  whole  of  the  pulse  as  it  propagates  and 
spreads. 

The  concept  of  the  parametric  array  was  introduced 
by  Westervelt  (1963)  and  led  to  considerably  increased 
interest  in  f i n i t e- amp  1 i tud e  acoustics.  A  parametric 
array  is  an  end-fire  array  formed  in  a  nonlinear  medium 
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by  Che  interaction  of  relatively  high-frequency  carrier 
beams  of  different  frequencies,  which  generate  a 
difference- frequency  signal.  As  the  amplitude  of  the 
carriers  is  reduced  by  absorption  and  spreading,  the  para¬ 
metric  array  has  a  built-in  amplitude  taper.  The  principal 
advantages  of  a  parametric  array  are  the  relatively  narrow 
beamwidth  of  the  d i f f e r e nee- f r e q ue n cy  signal,  which  is  com¬ 
parable  with  that  of  the  carriers,  the  virtual  absence  of 
side  lobes  in  the  difference-frequency  beam  in  many  cases, 
and  wide  relative  bandwidth.  The  principal  disavantage  is 
its  low  conversion  efficiency. 

An  asymptotic  far-field  value  for  the  difference 
frequency  level  arising  from  a  bifrequency  signal  has  been 
derived  by  Fenlon  (1972).  Let  the  acoustic  Reynold's  number 
be  defined  by 

r=3ek/a,  (1.10) 


where  £  is  the  Mach  number,  3  is  the  parameter  of 
nonlinearity,  k  is  the  wave  number  of  the  reference 
frequency,  in  this  case  the  difference  frequency,  and 
Ct=  6  £ 2  is  the  attenuation  coefficient  of  the  difference 
frequency.  The  downshift  ratio  Y0  is  the  ratio  of  the 
mean  carrier  frequency  to  the  difference  frequency.  A 
scaled  acoustic  Reynold's  number  may  be  defined  as 


r  =  r/2Yn 
1,2  0 


provided  that  f^~  f 


(1.11) 


,  i.  e.,  that  YQ  is  large,  and  that 
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the  carrier  amplitudes  are  equal;  then,  the  function 


P  Co)  =  1  I? (Ti , z/2)  -o/r 

y/i.2  TJTrTTTm 


(1.12) 


expresses  the  difference-frequency  level  as  a  function  of 
a  for  a>> 1 .  Here  ft  is  the  higher  of  the  two  carriers 
and  f2  is  the  lower,  with  0  the  scaled  range  as  defined 
previously.  Equation  (1.12)  may  be  multiplied  through  by 
exp(a/f)  so  as  to  cancel  the  effect  of  viscous  absorption, 


and  writ  ten  in 


p- (o/F) e" 


form 


1 1  (  F  i  ,  2  /  2  ) 
Iq(F 1,2/2) 


1  ( r  1 . 2  /  4 )  , 

i+(ri , 2/4) 2 


(1.13) 


where  the  latter  approximation,  involving  the  leading  terms 
of  the  Bessel  functions,  is  accurate  to  within  a  few 
percent.  The  right  side  of  this  equation  is  a  maximum  for 
F 1 ,2/4.  Since 

Ti , 2=  r/2yo,  (1.14) 

for  any  downshift  ratio,  there  is  a  value  of  T  that  will 
maximize  the  efficiency  of  the  parametric  conversion. 

The  value  of  T  for  which  maximum  efficiency  is  attained  is 


r 

ma  x 


(1.15) 


In  Fenlon  and  Mckendree' s  (1979)  solution  for  the 
propagation  of  a  bifrequency  signal  at  weak  finite 
amplitudes,  the  beam  shapes  of  the  carrier  waves  are 


approximated  by  Gaussian  functions,  and  it  is  assumed  that 
each  carrier  propagates  linearly.  Unlike  any  of  the  previous 
methods,  this  technique  includes  the  effects  of  diffraction 
(spreading  losses)  via  a  three-dimensional  form  of  Burger's 
equation  due  to  Zablotskaya  and  Kokhlov  (1969)  and  to 
Kuznetsov  (1971).  This  method  allows  the  difference- 
frequency  level  arising  from  interaction  of  the  carriers  to 
be  determined.  This  solution  applies  for  axisymmetric  waves 
propagating  in  a  viscous  fluid  and  can  be  evaluated  for  a 
wide  range  of  parameters.  However,  it  fails  at  large  source 
levels  due  to  its  neglect  of  finite-amplitude  losses  in  the 
carriers,  and  the  amount  of  the  error  cannot  be  determined 
as  a  function  of  the  source  level.  The  weak-finite-ampli¬ 
tude  solution  is  compared  with  the  results  of  the  numerical 
method  in  Chapter  4. 

Bakhvalov,  Zhileikin,  Zablotskaya,  and  Kokhlov  (1978, 
1979)  have  published  papers  to  date  on  finite  amplitude 
wave  propagation.  They  employ  a  direct  numerical  solution 
of  the  second-order  nonlinear  wave  equation,  but  their 
method  is  not  explained.  The  results  of  the  operational 
solution  are  compared  with  the  results  obtained  by 
these  researchers  in  Chapter  4. 
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Experiments  at  Finite  Amplitudes 

A  number  of  experiments  have  been  conducted  using 
sources  driven  at  finite  amplitudes.  Shooter,  Muir,  and 
Blackstock  (1974)  used  a  454  kHz  (kiloHertz)  source  in  fresh 
water.  Muir  and  Willette  (1972)  performed  experiments  with 
bifrequency  signals  having  carriers  in  the  vicinity  of 
450  kH2,  and  Eller  (1974)  used  carriers  near  1435  kHz. 

Each  of  these  experiments  has  been  simulated  by  the 
numerical  solution,  and  the  agreement  between  numerical 
and  experimental  results  is  discussed  in  Chapter  4. 
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Chapter  2 


DERIVATION 


This  chapter  presents  a  brief  derivation  of  the  second- 
order-nonlinear  wave  equation  in  the  form  in  which  it  will 
be  used  subsequently  in  this  thesis.  The  second  section 
reduces  the  equation  to  Burger's  form  and  shows  the  simpli¬ 
city  of  the  resulting  equation  in  scaled  coordinates.  The 
last  section  of  this  chapter  introduces  attenuation  and 
dispersion  into  Burger's  equation. 

The  Second-Order-Nonlinear 
Wave  Equation 


The  inviscid  form  of  the  nonlinear  wave  equation 
correct  to  second  order  as  given  by  Goldstein  (1960)  is 

(v2“  cT^*  “ 


-1(V^.V4>)+  |i  V2( 


(2.1) 


where  <P  is  the  velocity  potential.  The  first  term  within 
the  brackets  on  the  right  side  of  Equation  (2.1)  represents 
the  nonlinearity  introduced  by  convection  of  momentum 
through  the  element  of  fluid.  The  second  term  represents 
the  nonlinearity  of  the  equation  of  state  to  second  order. 
The  form  shown  is  for  an  ideal  gas;  in  a  fluid,  the  coeffi- 


in  the  form 


P ' 


P 

r  o  o 


(2.7) 


Scaled  Coordinates  and  Burger's  Equation 


In  nonlinear  equations,  a  system  of  scaled  coordinates 
is  used  that  emerges  naturally  from  normalization  of 
Equation  (2.7).  A  transformation  due  to  Kokhlov  (1961)  is 
applied  to  Equation  (2.7),  involving  the  dummy  range 
variable  £  and  substitution  of  the  retarded  time  t’: 


5  =  z,  T’=  T-z/cq, 


(2 .8a, b) 


where  is  the  direction  of  propagation.  If  the  definition 

9  2 


V2=  V2  +  — - 

I  9  z  ‘ 


(2.9) 


is  used,  then  the  Zab lo t ska ya-Kok hlo v  equation  results 
C 


.IlI  _  -2.  Vv  <  p>dT'= 


2 

_  9£J_. 

2  p  c  a  x  * 

r  o  o  ° 


32  -  T  ) 

whose  plane-wave  form  is  Burger's  equation: 

,  2 

9p 1  _  B  9  p  • 


(2.10) 


(2.11) 


9  z 


2  p  c  0  9  i  • 

o  o 


To  obtain  this  result,  it  is  assumed  that  the  rate  of 


change  of  the  waveform  with  respect  to  distance  is  very 
much  less  than  its  rate  of  change  with  respect  to  time. 


In  finite-amplitude  propagation  in  fluids,  this  assumption 
is  justified  as  the  cumulative  effects  of  nonlinearity 
occur  over  many  wavelengths  of  the  signal. 


Now  the  retarded  time,  range,  and  pressure  are 
normalized  with  respect  to  the  subscript-star  references: 

T  =  lo^t'  =  w^(t-z/co)  =  aj^t-k^z,  (2.12a) 

a  =  Be*k*z  (2.12b) 

and 

P  =  P’/P*  (2.12c) 

where  ,  k*,£*  ,  and  p^  are  the  frequency,  wavenumber, 

Mach  number,  and  peak  pressure  of  an  arbitrary  reference 
wave.  Substitution  of  the  above  into  Burger's  equation 
leads  to  the  equation  in  scaled  coordinates: 

_3P  =  1  3P2  (2.13) 

3a  2  3x 


Attenuation  and  Dispersion 

The  presence  of  viscous  forces  and  of  dispersivity  in 
the  medium  introduces  additional  terms  into  the  forces 
acting  on  each  element  of  fluid.  These  effects  are  linear 
and  thus  do  not  affect  the  nonlinear  term  as  represented 
in  Burger's  equation. 

Attenuation  of  sound  in  a  fluid  is  produced  by  viscous 


losses  and  by  thermal  relaxation 


The  t hermoviscous 


attenuation  parameter  of  a  fluid,  a,  is  defined  by  Landau 
and  Lifshitz  (1959)  as 


r  =  2n+n'  ,  k(y-i)/cp 

2p  c  2  2p  c  2 

O  O  K  O  O 


(2.14) 


in  which  p  and  p'are  the  shear  and  dilatational  coefficients 
of  viscosity,  p^  and  are  the  equilibrium  pressure  and 
density,  K  is  the  thermal  conductivity,  Y  is  the  ratio  of 
specific  heats,  and  Cp  is  the  specific  heat  at  constant 
pressure.  The  first  part  of  Equation  (2.14)  represents 
viscous  absorption,  and  the  second  represents  thermal 
relaxation.  In  all  fluids  except  liquid  metals,  the  thermal 
relaxation  term  is  several  orders  of  magnitude  smaller  than 
the  viscous  absorption  term. 

The  thermoviscous  attenuation  coefficient  of  a  fluid  at 
the  frequency  f  is  defined: 


a=  6  f ' 


(2.15) 


and  in  progressive  plane-wave  propagation,  the  effect  of 
absorption  is  an  exponential  decrease  in  amplitude  as  a 
function  of  the  distance  traveled: 


p  (z)  =  p (o  )  e 


-az 


(2.16) 


which  applies  for  the  value  of  Ct  corresponding  to  each 
frequency  component. 

In  fluids  which  contain  or  are  made  up  of  polyatomic 
molecules,  different  molecular  configurations  may  have 
different  internal  energies.  It  is  possible  for  acoustic 
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energy  to  be  stored  in  a  rearrangement  of  the  molecular 
structure  and  subsequently  restored  to  the  acoustic  signal. 
Each  internal  relaxation  process  is  characterized  by  its 
relaxation  time  T  .  A  monorelaxing  fluid  exhibits  the 
property  of  dispersion — that  is,  a  linear  dependence  of  the 
velocity  of  propagation  on  frequency — at  frequencies  far 
removed  from  its  relaxation  frequency.  The  dispersivity 
of  a  fluid  is  defined 

®  =  -  c20  )/c>  •  (2.17) 

The  complete  derivation  of  the  second-order-nonlinear 
wave  equation  including  the rmovis cous  absorption  and 
relaxation  terms  is  outside  the  scope  of  this  thesis. 

The  interested  reader  is  referred  to  the  derivation  by 
Lamb  (1959)  for  the  dispersion  term  and  the  wave 
equation  derivation  by  Kuznetsov  (1971). 

For  a  monorelaxing  fluid  having  relaxation  time  and 
attenuation  parameter  a,  the  acoustic  Reynold's  number  may 
be  wr i t  ten : 


T  =  Se*k^/a,  (2.18) 

and  the  scaled  relaxation  and  dispersion  coefficients  A 
and  d  may  be  written 

=  28e^/m 


(2.19) 
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and 


d  =  2ee*/m(w*TR)2, 


(2.20) 


each  evaluated  at  the  scaling  frequency  m  and  Mach  number 

* 

r  .  In  scaled  coordinates,  the  s e co nd -o r d e r- non  1 in e a r  wave 
equation  including  the  effects  of  relaxation,  thermoviscous 
losses,  and  diffraction  takes  the  form  of  an  equation  first 
derived  in  slightly  different  form  by  Kokhlov  (1961).  In  a 
non-relaxing  fluid,  this  equation  reduces  to  Kuznetsov's 
equation  (1971)  and,  in  a  non-relaxing  and  lossless  fluid, 
it  reduces  to  the  Zab lo t s kay a-Kokh lo v  equation  (1969).  The 
second-order-nonlinear  wave  equation  in  scaled  coordinates 


including  relaxation,  spreading,  and  losses,  takes  the  form 


(1+“*Tr  fr> 


I  a2p 
r  3t2' 


1  9  p2 

2  3  t 


a)* TR  3  2  P 

A  3  t  2  (2.21) 


where  the  transverse  Laplacian  operator  is  primed  to 

indicate  that  it  is  in  scaled  coordinates.  If  <<1, 

R 

then  Equation  (2.21)  takes  a  simpler  form  including 
dispersion  rather  than  relaxation: 
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i  a 3p  l  a2p 

d  at3  ~  r  aT7 


1  3P2 

2  3x 


(2.22) 


In  Equation  (2.22),  the  nonlinearity  is  represented  in 
the  same  form  as  Burger's  equation.  The  second  term 
represents  dispersion,  the  third  represents  thermoviscous 
absorption,  and  the  fourth  represents  diffraction.  The 


plane-wave  form  of  Equation  (2.22)  is 
3j?  _  i  3  3P  1  32P  =  i  ail 

3  a  d  3x 3  ”  T  FP-  2  3  T  »  (2.23) 

which  is  the  limit  as  tends  to  infinity,  suppressing 
the  diffraction  term.  In  a  lossless  fluid,  the  acoustic 
Reynold's  number  tends  to  infinity,  removing  the  absorption 
term.  In  a  non-relaxing  fluid  m  is  zero,  so  d  becomes 
infinite  and  removes  the  dispersion  term. 

Equations  (2.21),  (2.22),  and  (2.23)  were  derived  using 
the  progressive-wave  transformation  due  to  Kokhlov,  as  was 
used  in  the  scaling  of  Equations  (2.10)  and  (2.11);  thus, 
the  same  restrictions  apply.  In  addition,  a  paraxial 
assumption  requires  that  the  effects  of  attenuation  and 
dispersion  are  significant  only  along  the  direction  of 
propagation;  the  transverse  components  of  attenuation 
and  dispersion  are  ignored. 
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Chapte  r  3 


NUMERICAL  SOLUTION 


This  chapter  presents  each  feature  of  the  numerical 
solution  separately.  The  second-order-nonlinear  wave 
equation.  Equation  (2.22),  is  solved  in  parts:  first  the 
lossless  Burger's  equation;  then  that  part  representing  loss 
and  dispersion;  and  last,  the  diffraction  term.  It  is 
assumed  that  each  of  these  effects  operates  independently  of 
the  others  over  small  distances,  so  that  it  is  possible  to 
consider  them  one  at  a  time.  That  is,  a  sampled  wavform  is 
propagated  over  an  incremental  range  according  to  the 
lossless  Burger's  equation.  It  is  then  Fourier  transformed, 
and  each  frequency  component  is  subjected  to  that  attenu¬ 
ation  and  dispersion  which  it  would  have  suffered  in  propa¬ 
gating  over  the  same  distance  in  a  linear  manner.  For 
axisymmetric  three-dimensional  propagation,  the  procedure 
outlined  above  is  applied  to  each  of  a  set  of  sampled  wave¬ 
forms  at  various  radial  distances;  then  the  resulting 
signal  is  subjected  to  linear  diffraction  over  the 
incremental  distance. 

The  first  section  of  this  chapter  discusses  the 
numerical  solution  of  the  lossless  Burger's  equation.  The 
next  section  presents  the  method  by  which  attenuation  and 


dispersion  are  introduced  in  the  frequency  domain,  and 
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the  third  section  explains  the  numerical  diffraction 
method.  The  last  section  of  this  chapter  discusses 
features  of  the  computer  programs  which  implement  the 
numerical  solution,  and  reviews  the  reasons  for  certain 
choices  made  in  the  methods  of  implementation. 

Lossless  Plane-Wave  Solution 


The  lossless  plane-wave  solution  is  obtained  by  direct 
application  of  a  method  due  to  Bellman  et  al.  (1965).  In 
this  method  the  points  on  the  waveform  advance  at  a  velocity 
that  is  proportional  to  their  amplitudes.  Beginning  with 
the  lossless  Burger's  equation, 


_  P  3iL  =  0  ,  (3.1) 

3a  3t 

the  derivative  with  respect  to  a  may  be  approximated  as 


|1  n,  i 
3  a  =  a 


P (a , x )  -  P  (0,x) 


(3.2) 


where  a  is  used  in  the  sense  of  an  incremental  range. 
Equations  (3.1)  and  (3.2)  become,  by  rearranging  terms, 


r\  P 

P  (a, t)  =  P(0,t)  +  oP(o,t)3~ 


(3.3) 


If  the  quantity  oP(o,t)  is  small,  then  it  is  reasonable  to 
relate  Equation  (3.3)  to  a  Taylor's  series  expansion.  If  p 
is  slowly  varying  then,  correct  to  second-order  terras, 

P(0,t)  =  P 


— 

C 


0 , T+OP (0 , T ) 


(3.4) 
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This  formulation  is  not  subject  to  failure  by  prediction 
of  multivalued  pressures,  as  each  point  of  the  waveform  at 
the  advanced  range  is  related  to  a  definite  point  at  the 
initial  range.  Other  formulations  involving  explicit 
numerical  differentiation  in  a  form  reminiscent  of 
Equation  (3.3)  have  been  found  to  become  unstable  in  the 
vicinity  of  a  discontinuity. 

The  numerical  implementation  may  be  enhanced  by  inclu¬ 
ding  an  estimate  of  P(cr,T)in  the  expression  for  the  incre¬ 
mental  time.  Using  Equation  (3.4),  an  estimate  of  P(o,t)  , 
denoted  P  (a,t),  may  be  defined  as 


Pg (a, t) 


T+aP (0 , t) 


(3.5) 


and  this  estimate  may  be  substituted  into  Equation  (3.4) 
to  give 


P (a, t) 


0, 


T+aP  (o,t) 
e 


) 


(3.6) 


Equations  (3.4)  and  (3.6)  are  implemented  in  the  fol¬ 
lowing  manner.  A  number  of  samples  of  the  pressure  waveform 
are  taken  at  equal  intervals  covering  the  data  window. 

These  samples  represent  P(0,T).  A  second  array  of  samples 
is  prepared,  each  of  whose  values  is  interpolated  from  among 
the  P(0,T)at  the  time  specified  by  Equation  (3.4).  These 
are  the  values  of  the  waveform  at  the  range  o.  The  values 
of  the  waveform  at  the  new  range  may  be  used  as  the  origin 
for  a  new  range  step,  and  the  procedure  may  be  repeated 
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until  the  desired  range  is  reached.  Both  equations  are 
based  on  a  Taylor's  series  and  are  valid  for  o  less  than  1. 
For  implementation  of  Equation  (3.6),  the  first  estimate 
is  used  to  obtain  a  new  estimate  of  the  advanced  waveform 
at  each  range  step. 

Attenuation  and  Dispersion  in  the 
Frequency  Domain 


If  the  normalized  pressure  waveform  P  is  Fourier 
transformed,  its  spectrum  will  consist  of  the  coefficients 
P  at  frequencies  which  are  related  by  integers. 

GJ 

Equation  (2.22)  expressed  in  the  frequency  domain  takes 
the  form 


9  Pu)  ,  ^ 
9a  C  d 


■)Pw  =  2^^*  * 


(3.7) 


Without  loss  of  generality,  it  may  be  assumed  that  u)  is 
an  integer  multiple  of  ,  the  frequency  at  which  the 
scaling  of  Equations  (2.21)  et  seq.  was  performed.  Then, 
the  equation  in  the  frequency  domain  becomes 


2_En  + 

9  a 


AB.p  *  p 
2  n  n 


(3.8) 


with  the  n-th  Fourier  harmonic  component. 

For  a  mono-relaxing  fluid  having  relaxation  time  , 
additional  terms  in  the  scaled  f  requency  nu^  are  introduced 
by  the  derivatives  in  Equation  (2.21).  The  frequency-domain 
equation  then  takes  the  form 


2  GJ  T  r»  /  * 

(D_  ,,  *_  R r/A  + 

lT  1+  (ruo^T  R  )  i 


1+(nw*TR) 


=  An  P  *P  . 
— 2  n  n 


(3.9) 


Since  the  method  of  the  first  section  provides  a  satis¬ 
factory  solution  of  the  lossless  equation,  it  may  be  used  to 
propagate  the  waveform  to  an  incrementally  advanced  range. 
The  waveform  is  then  Fourier  transformed,  and  each  of  its 
spectral  components  is  processed  according  to  the  linear 


equation  for  attenuation  and  dispersion: 

9Pn  (  ini  +  ni  )p 

3c  1  d  T  ;  n  * 


(3.10) 


for  which  the  formal  solution  is 


1  i  2 

P  (ct,t)  =  P  (0,T)e  v  d  T  ;  • 
n  n 


(3.11) 


Note  that  the  coefficient  of  ?n  in  Equation  (3.10)  is 
transferred  intact  into  the  exponential  in  the  solution. 
Thus,  the  effects  of  an  arbitrary  attenuation  and  dispersion 
as  a  function  of  frequency  may  be  introduced  by  multiplying 
the  spectrum  by  a  complex  function  of  frequency  whose 
form  may  be  determined  from  the  governing  wave  equation 
by  inspection. 


Diffraction  by  Diffusion  in  Frequency 


The  nonlinear  wave  equations,  Equation  (2.21)  et  seq., 


contain  diffraction  terms. 


If  the  other  terms  are  taken 


into  account  by  the  methods  of  the  preceding  sections,  then 


the  remaining  equation  takes  the  form  of  a  diffusion 
equation  in  the  frequency  domain: 


+  _i_ 

3  a  no 

o 


0. 


(3.12) 


The  effects  of  diffraction  on  the  propagation  of  finite 
amplitude  waves  in  three  dimensions  is  introduced  by  a 
finite-difference  numerical  solution  of  Equation  (3.12); 
this  method  is  termed  "diffusion  in  frequency". 

The  incremental  range  O  may  be  chosen  small  enough  to 
render  negligible  all  terms  higher  than  the  first  order. 
Under  these  conditions,  Equation  (3.12)  may  be  converted  to 


pn(a)  *  Pn(o) 


^  \7 2 

na0  ^Vo)’ 


(3.13) 


which  is  easily  solved  by  f ini t e-di f f e r enc e  methods. 

A  finite-difference  numerical  technique  must  be  used  as  the 
functional  form  of  pn  is  not  known. 


Comments  on  the  Numerical  Solution 


The  method  of  numerical  solution  of  the  second-order- 
nonlinear  wave  equation  described  in  this  chapter  is  imple¬ 
mented  by  a  group  of  computer  programs  written  in  ANSI 
standard  FORTRAN  IV  and  operating  on  an  IBM  370/3033 
running  under  OS/MVT.  On  this  system,  data  sample  arrays  of 
almost  unlimited  size  may  be  accommodated.  The  data  arrays 
in  the  program  are  restricted  to  236  by  40  points  so  that 
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Che  program  will  fit  into  a  smaller  computer,  also  used 
during  this  research.  The  processing  time  for  simulations 
which  exercise  the  program  most  severely  is  on  the  order  of 
300  to  1000  seconds.  Plane  wave  propagation  may  be  modeled 
with  up  to  2048  waveform  samples.  Such  simulations  require 
comparable  lengths  of  time  for  execution. 

The  evaluation  of  Equation  (3.4)  or  (3.6)  requires 
interpolation  between  samples  of  the  waveform.  Several 
methods  have  been  tried,  including  spline  fitting,  second- 
through  fifth-order  fitted  polynomials,  and  third-  through 
seventh-order  Lagrange  interpolation.  Each  of  these 
advanced  methods  fails  in  the  vicinity  of  the  discontinuity. 
Accordingly,  a  linear  fit  interpolator  is  used.  The  advanced 
time  for  each  new  waveform  sample  falls  between  two  samples 
of  the  old  waveform,  and  the  interpolated  value  is  assumed 
to  lie  on  the  straight  line  joining  the  old  samples.  This 
procedure  has  been  found  to  be  sufficiently  accurate  by 
comparison  with  known  analytic  solutions  in  their  domains 
of  applicability.  The  method  is  biased,  however,  in  the 
following  manner.  A  sinusoidal  signal  is  everywhere  con¬ 
cave  with  respect  to  the  ordinate  axis.  A  linear  inter¬ 
polator  will  therefore  tend  to  underestimate  its  values, 
as  the  sinusoid  lies  further  from  the  ordinate  axis  between 
nearby  samples,  than  does  the  straight  line  joining  them. 

This  property  of  the  linear  interpolator  introduces  an 
error  or  bias  into  the  numerical  results.  Other  waveforms 


of  practical  interest  have  segments  which  are  convex  with 
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respect  to  the  ordinate  axis  and  are  there  overestimated. 


The  amount  of  the  error  may  be  reduced  by  taking  more  wave¬ 
form  samples.  The  accuracy  of  the  interpolation  is  also 
affected  by  the  number  of  steps  per  unit  range.  If  the 
step  size  is  small,  the  advanced  times  at  which  the  inter¬ 
polated  waveform  values  are  desired  will  be  close  to  the 
times  of  the  existing  waveform  samples,  and  the  error  of  the 
linear  interpolation  will  be  less  than  if  the  advanced  times 
lay  further  from  the  existing  samples. 

It  is  instructive  to  consider  the  sample  number  con¬ 
straints  in  the  frequency  domain.  A  sequence  of  N  samples 
spanning  a  time  T  can  be  Fourier  transformed  into  the  ampli¬ 
tudes  of  frequency  components  which  are  the  harmonics 
of  1/T.  it  has  been  found  that  acceptable  results  are 
obtained  in  the  numerical  solution  if  at  least  ten  harmonics 
of  the  highest  significant  frequency  component  at  the  source 
can  be  accommodated  in  the  spectrum,.  A  sequence  of  length  N 
yields  a  spectrum  which  coversN/2+1  harmonics  of  1/T. 

Thus  N/2+1  must  be  chosen  greater  than  or  equal  to  ten 
times  the  relative  frequency  of  the  highest  frequency  source 
component.  For  a  pure  tone  input,  the  highest  significant 
source  frequency  is  the  input  tone.  For  a  bifrequency 
signal,  the  sum  frequency  must  have  about  ten  harmonics  to 
ensure  reasonably  good  results  from  the  nonlinear 
operator. 

Certain  kinds  of  signals,  for  example,  trains  of  soli- 
tons ,  will  propagate  in  an  appropriate  medium  without  loss 
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of  energy.  For  this  reason,  the  numerical  solution  as 
implemented  in  the  computer  programs  includes  an  optional 
feature  which  ensures  conservation  of  energy  in  the  non¬ 
linear  operator.  This  option  must  be  used  with  discretion, 
as  many  kinds  of  signals  are  subject  to  finite-amplitude 
losses. 

For  R>>1,  that  is,  at  ranges  much  greater  than  the 
Rayleigh  distance,  spherical  spreading  may  be  assumed.  The 
diffraction  by  diffusion-in-frequency  operator,  as  imple¬ 
mented  in  its  subroutine,  allows  the  user  to  select  a  value 
of  R  beyond  which  spherical  spreading  is  assumed.  Since  the 
Rayleigh  distance  differs  for  each  harmonic,  the  values  of  ° 
at  which  spherical  spreading  may  be  imposed  differ  in  a 
corresponding  way,  in  direct  proportion  to  their  Fourier 
harmonic  numbers. 

The  effects  of  diffraction  may  be  taken  into  account 
by  any  of  several  different  techniques.  It  is  possible  to 
represent  the  beam  function  at  each  frequency  as  a  sum  of 
eigenmodes  of  the  geometry  in  use;  then,  the  diffraction  may 
be  interpreted  as  a  phase  modulation  of  the 
spatial-frequency  spectrum.  Attempts  to  implement  this 
method  by  numerical  means  were  usable  but  subject  to 
breakdown  in  the  far  field.  The  finite  number  of  samples 
in  the  radial  direction  implies  a  spatial  periodicity  ol 
the  source  function,  and  interference  between  the  image 
sources  becomes  evident  in  the  far  field. 


Chapter  4 


RESULTS 


This  chapter  presents  results  of  the  application  of  the 
numerical  solution  to  a  number  of  problems.  The  first 
section  gives  examples  of  the  independent  operation  of 
certain  parts  of  the  numerical  solution  to  demonstrate  its 
degree  of  stability  and  accuracy.  The  second  section  deals 
with  finite  amplitude  propagation  of  plane  waves,  and 
discusses  the  coupled  effects  of  nonlinearity,  attenuation, 
and  dispersion  as  they  are  introduced  in  various  relative 
strengths.  The  third  section  shows  the  results  of  the 
weak  finite  amplitude  solution  and  those  of  the  numerical 
solution  for  equivalent  conditions,  including  diffraction 
and  spreading.  The  fourth  section  presents  the  results  of 
certain  experiments  at  finite  amplitudes  and  their 
numerical  simulations.  The  last  section  presents  numerical 
simulations  of  various  types  of  modulated  signals  in 
axisymmetric  propagation  in  one  and  in  three  dimensions  at 
finite  amplitudes. 


Numerical  Solution  of  the  Lossless  Burger's 
Equation  and  of  Linear  Diffraction 


The  number  of  waveform  samples  and  the  number  of  steps 


per  unit  of  range  may  be  chosen  at  the  discretion  of  the 


user  of  the  numerical  solution.  Table  4.1  presents  the 
amplitudes  of  the  first  three  harmonics  of  a  pure-tone 
signal  at  two  ranges  as  a  function  of  the  number  of 
waveform  samples,  with  10  steps  per  unit  range.  Table  4.2 
presents  the  amplitudes  of  the  first  three  harmonics  of  a 
pure-tone  signal  at  two  ranges  as  a  function  of  the  number 
of  steps  per  unit  range  for  64  waveform  sample  points. 

For  each  table,  the  analytic  value  according  to  the  Fay 
solution  and  the  error  as  a  percentage  of  this  value 
is  shown.  These  results  were  derived  from  the  use  of 
Equation  (3.4)  after  the  method  due  to  Bellman  (1965). 

If  the  revised  expression  given  by  Equation  (3.6)  is 
used,  then  better  results  may  be  obtained  with  a  smaller 
number  of  calculations.  Table  4.3  shows  the  amplitudes  of 
the  first  three  harmonics  of  a  pure-tone  signal  at  various 
ranges,  as  deduced  from  the  matched-asymptotic  solution  due 
to  Blackstock  (1966)  and  from  the  numerical  implementation 
of  Equation  (3.6).  In  this  case  128  waveform  samples  were 
used,  with  10  steps  per  unit  o  to  0=1  and  5  steps  per  unit 
O beyond  0=1.  Errors  are  shown  between  the  best  numerical 
results  and  the  matched-asymptotic  solution  values  as  a 
percentage  of  the  latter. 

Figure  1  shows  the  level  on  axis  of  a  Gaussian  beam, 
exp(-e2  )  at  R  =0 ,  and  Figure  2  shows  the  beam  widths  to 
the  -3,  -6,  and  -10  dB  points,  as  measured  in  the  numer¬ 
ical  solution  as  a  function  of  range.  The  discrete  points 


are  derived  from  the  implementation  of  the  d i f f u s  i  on- i n- f r e- 
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Table  4.1.  Harmonic  Levels  at  Two  Ranges  for  an  Initially 
Pure-Tone  Signal  Propagating  in  a  Lossless 
Nonlinear  Medium,  For  Various  Numbers  of 
Waveform  Samples  and  10  Steps  Per  Unit  Range, 
Derived  from  the  Numerical  Solution  of 
Equation  (3.4). 


NUMBER  OF 


WAVEFORM 

B  (5) 

B  (5) 

B  (5) 

B  (10) 

B  (10) 

B  (10) 

SAMPLES 

1 

2 

3 

1 

2 

3 

8 

.  288 

.  108 

.040 

.  158 

.058 

.020 

16 

.310 

.  142 

.083 

.  170 

.076 

.042 

32 

.317 

.153 

.096 

.  174 

.082 

.050 

64 

.320 

.157 

.102 

.  175 

.084 

.053 

128 

.321 

.  158 

.  104 

.  176 

.  086 

.055 

256 

.321 

.159 

.  104 

.  176 

.  086 

.056 

512 

.321 

.  159 

.  104 

.  176 

.086 

.056 

1024 

.321 

.  159 

.  104 

.  176 

.  086 

.056 

analytic 

.333 

.  167 

•  111 

.  182 

.091 

.061 

error,  % 

-3.6 

-4.8 

-6.3 

-3.3 

-5 . 5 

-8.2 

Table  4.2.  Harmonic  Levels  at  Two  Ranges  for  an  Initially 
Pure-Tone  Signal  Propagating  in  a  Lossless 
Nonlinear  Medium,  For  Various  Number  of  Steps 
Per  Unit  Range  and  64  Waveform  Samples, 

Derived  from  the  Numerical  Solution  of 
Eq  ua t ion  (3.4). 

NUMBER  OF 


RANGE 

B  (5) 

B  (5) 

B  (5) 

B  (10) 

B  (10) 

B  (10) 

STEPS 

1 

2 

3 

1 

2 

3 

5 

.311 

.153 

.099 

.  172 

.083 

.052 

10 

.320 

.157 

.  102 

.175 

.084 

.053 

20 

.324 

.  159 

.  103 

.  177 

.086 

.054 

40 

.326 

.  160 

.  103 

.  178 

.086 

.054 

analytic 

.333 

.  167 

.111 

.  182 

.091 

.061 

error ,  % 

-2.  1 

-4 . 2 

-6.3 

-2.2 

1 

Ln 

V-n  1 

-11.5 
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Table  4.3.  Levels  of  the  First  Three  Harmonics  of  a  Pure 
Tone  Signal  in  Lossless  Nonlinear  Plane-Wave 
Propagation  as  Derived  from  the  Matched-Asymp¬ 
totic  Solution  due  to  Blackstock  (1966)  and 
from  the  Numerical  Solution  of  Equation  (3.6), 
with  Errors  Expressed  as  a  Percentage  of  the 


Former . 

MATCHED- 

NUMERICAL 

ERROR 

SIGMA 

ASYMPTOTIC 

EQUATION  (3.6) 

PERCEN 

B  (a) 

B  ( a ) 

1 

1 

0.1 

0.9988 

0.9986 

-0 . 02 

0.2 

0.9950 

0.9948 

-0.02 

0.5 

0.9691 

0.9685 

-0.06 

1  .0 

0.8801 

0.8820 

0.22 

2.0 

0.6484 

0.6518 

0.52 

5.0 

0.3333 

0.3349 

0.48 

10.0 

0.1818 

0. 1816 

-0.11 

B  (  a) 

B  (a) 

2 

2 

0.1 

0.0498 

0.0496 

-0.40 

0.2 

0.0987 

0.0981 

-0.61 

0.5 

0.2298 

0.2284 

-0.61 

1 .0 

0.3528 

0.3482 

-1.3 

2.0 

0.3134 

0.3113 

-0 .67 

5.0 

0.  1667 

0.1679 

0.72 

10.0 

0.0909 

0.0932 

2.5 

B  (a) 

B  (a) 

3 

3 

0.1 

0.0037 

0.0037 

0.0 

0.2 

0.0147 

0 -0144 

2.0 

0.5 

0.0802 

0.0803 

0.12 

1  .0 

0 . 2060 

0. 2023 

-1.8 

2.0 

0.2070 

0.2047 

-1  .  1 

5.0 

0.1111 

0.1117 

0.54 

10.0 

0.0606 

0.0584 

-3 . 6 

i 


Source 


id  th 
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quency  technique  which  is  used  in  the  numerical  solution. 

The  known  analytic  values  for  the  diffraction  of  a  Gaussian 
beam  are  shwon  in  these  figures  as  solid  curves.  The  accu¬ 
racy  with  which  the  numerical  solution  matches  the  analytic 
values  indicates  that  the  numerical  solution  of  the  dif¬ 
fraction  term  by  diffusion-in-frequency  is  satisfactory. 

The  remaining  difficulty  with  the  diffusion-in-frequency 
technique  will  be  discussed  in  Chapter  5. 

Finite-Amplitude  Propagation 
of  Plane  Waves 

The  finite-amplitude  propagation  of  plane  waves 
in  a  thermoviscous  mono-relaxing  fluid  is  governed  by  three 
effects:  nonlinearity,  attenuation,  and  dispersion.  This 
section  demonstrates  the  operation  of  the  numerical  solution 
of  Equation  (2.22)  in  each  of  these  effects,  taken  sepa¬ 
rately  and  in  combination. 

The  nonlinear  operator  may  be  tested  by  comparison  with 
Blackstock's  matched  asymptotic  solution  to  the  lossless 
Burger's  equation.  Figure  3  shows  the  harmonic  amplitudes 
which  are  predicted  by  Blackstock's  formula;  Figure  4  shows 
the  same  function  as  predicted  by  the  plane-wave  lossless 
numerical  solution.  The  difference  between  the  matched- 
asymptotic  and  the  numerical  solution  is  nowhere  greater 
than  a  few  percent;  this  indicates  the  accuracy  of  the 
numerical  solution  in  lossless  plane-wave  propagation. 


.....  -  im  r  ‘  lit  I 


Amplitude  Relative  to  Peak  at  the  Source 
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scaled  Range ,  a 


Figure  3.  Lossless  Plane-Wave  Harmonic  Coefficients 
from  the  Matched  Asymptotic  Expression 
Derived  by  Blackstock. 


Amplitude  Relative  to  Peak  at  the  Source 


0  . 


0  . 


0. 


0. 


Figure  4  . 


Lossless  Plane- 
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Viscous  losses  introduce  an  exponential  decay  of  signal 
level  with  range.  WithP^  the  n-th  harmonic  amplitude, 

0  the  range,  and  T  the  attenuation  parameter, 


P 

n 


(a) 


P 

n 


(  0  )  e 


-n  2a/  P 

> 


(4.1) 


is  an  obvious  consequence  of  Equation  (3.10)  with  respect 
to  viscous  losses. 

Figure  3  shows  the  waveform  of  a  monofrequency  signal 
at  a  =1,  2,  and  3,  as  it  propagates  in  a  lossless  medium. 

The  waveform  steepens  until  a  shock  is  formed,  and  the  shock 
wave  decays  due  to  finite  amplitude  loss  without  changing 
its  shape.  In  Figure  6,  attenuation  has  been  introduced 
with  an  acoustic  Reynold's  number  of  10.  The  amplitude  of 
the  wave  is  reduced  by  finite-amplitude  losses,  and  in 
addition  each  harmonic  is  attenuated  by  the rmov i s co us 
effects.  Figure  7  shows  the  waveforms  with  very  high 
attenuation,  corresponding  to  r  =3.  At  O  =3,  the  waveform  is 
nearly  restored  to  a  sinusoid.  Table  4.4  lists  the  levels 
of  the  first  three  harmonics  of  each  of  the  waveforms  at 
each  range.  These  results  show  the  profound  effect  on  the 
waveform  and  spectrum  which  is  produced  by  viscous  absorp¬ 
tion.  As  the  rate  of  thermoviscous  absorption  is  usually 
proportional  to  the  square  of  the  frequency,  the  effect  on 
higher-frequency  components  is  even  more  marked. 

Modulated  waveforms  are  also  subject  to  shock 
formation.  Figures  8  through  10  show  a  bifrequency  signal 
propagating  in  a  lossless  medium.  The  large  amplitude 


i 

i 


d 


Relative  to  Peak  at  the  Source 


Table  4.4.  Amplitudes  of  the  First  Three  Harmonics  of  a 

Monofrequency  Signal  in  Plane  Wave  Propagation, 
with  Selected  Amounts  of  Absorption,  at 
Selected  Ranges. 


sigma 

B  (a) 

B  (a) 

B  (a) 

1 

2 

3 

(lossless) 

0.5 

.  96b2 

.  2292 

.  0786 

1.0 

.8760 

.  3525 

.  2045 

2.0 

.  6427 

.  3097 

.2038 

3.0 

.4906 

.  242  1 

.  1605 

5.0 

.  3300 

.  1648 

.  1097 

10.0 

.  1805 

.0907 

.  0598 

(  GAMMA= 10 ) 

0.5 

.9219 

.  1988 

.0600 

1  .0 

.8133 

.  2837 

.  1362 

2.0 

.6042 

.  2635 

.1514 

3.0 

.4619 

.  2042 

.1152 

5.0 

.  3055 

.  1248 

.0621 

(GAMMA=3 ) 

0.5 

.8253 

.  1423 

.0320 

1.0 

.6684 

.  1585 

.0462 

2.0 

.4439 

.  1055 

.0281 

3.0 

.  3043 

.0586 

.0120 

5.0 

.  1508 

.0162 

.0018 

Amplitude  Relative 


Amplitude  Relative  to  Peak  at  the  Source 


Figure  9 


Amplitude  Relative  to  Peak  at  the  Source 


cycles  go  into  shock  quickly  and  are  reduced  by  finite 
amplitude  effects  more  rapidly  than  the  small  amplitude 
cycles.  The  result  is  distortion  of  the  envelope  function. 
Figure  11  shows  an  envelope  function  constructed  by 
connecting  the  peaks  of  the  cycles  of  the  modulated 
waveform.  As  the  wave  propagates,  the  distortion  of  the 
envelope  becomes  more  evident.  Figure  12  shows  modulation 
envelopes  for  a  larger  downshift  ratio. 

Bifrequency  signals  propagating  in  a  nonlinear  medium 
give  rise  to  combination  frequencies  due  to  the  self¬ 
convolution  of  the  signal.  If  the  two  carriers  are  denoted 
fi  and  f2  ,  the  lowest  frequency  in  the  nonlinearly  generated 
spectrum  is  the  difference  frequency  f=fj-f2  •  Figure  13 
shows  the  difference-frequency  level  as  a  function  of  O 
generated  by  a  bifrequency  signal  for  each  of  several 
values  of  T.  After  all  f i ni t e-amp  1 i t ud e  effects  have  ended, 
the  plane-wave  difference-frequency  signal  will  continue  to 
propagate  subject  only  to  t he rmovi s cous  losses.  Figure  14 
shows  the  d if f e r e nee- f r e q ue ncy  level  multiplied  by 
exp(t?/  T)  ,  Each  of  the  curves  is  tending  towards  a  limiting 
value.  Figure  13  shows  the  far-field  difference-frequency 
level  times  exp(a/D  as  a  function  of  F  for  three  downshift 
ratios.  The  discrete  points  are  derived  from  the  numerical 
solution,  and  the  solid  curves  are  from  the  asymptotic 
far-field  solution  due  to  Fenlon  (1972).  The  results 
derived  from  the  numerical  solution  are  very  similar  to  the 
exact  analytic  results  as  either  the  acoustic  Reynold's 


Retarded  Time,  T 


Modulation  Envelopes  of  a  Bifrequency 
Signal  Having  a  Frequency  Downshift  Ratio 
of  5.5,  at  Selected  Ranges. 


Retarded  Time,  T 


Figure  12 


Modulation  Envelopes  of  a  Bifrequency 
Signal  Having  a  Frequency  Downshift  Ratio 
of  10.5,  at  Selected  Ranges. 


Acoustic  Reynold's  Number,  i 


Far-Field  Difference-Frequency  Level 
Times  exp(o/r),  as  a  Function  of  F, 
for  Selected  Frequency  Downshift 
Ratios. 
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number  or  the  frequency  downshift  ratio,  or  both,  are 
varied.  That  correct  results  are  obtained  as  either  is 
varied  implies  that  both  the  scaling  of  the  nonlinear  oper¬ 
ator  with  respect  to  frequency,  and  the  effects  of  thermo- 
viscous  absorption,  are  correctly  modeled. 

The  presence  of  nonlinearity  in  a  lossless  dispersive 
medium  gives  rise  to  the  formation  of  solitons.  Figure  16 
shows  the  numerical  simulation  of  the  interaction  of  two 
solitons  as  they  propagate  past  one  another.  The  solitons 
interact  and  pass  through  each  other  without  change.  By 
contrast,  a  sinusoidal  input  signal  in  a  nonlinear  and 
dispersive  medium  will  be  resolved  into  a  set  of  solitons. 
Each  of  these  will  propagate  at  its  own  speed.  If  a  con¬ 
tinuous  wave  signal  is  thus  resolved  into  a  soliton  train, 
it  is  possible  for  the  solitons  to  coalesce  at  later  times 
and  to  regenerate  the  initial  waveform.  This  behavior  has 
been  deduced  on  analytic  grounds  by  Zabusky  and  Kruskal 
(1965),  and  is  illustrated  in  the  following  two  figures. 

Figure  17  shows  the  history  of  an  initially  sinusoidal 
signal  as  it  propagates  through  a  lossless,  nonlinear,  and 
dispersive  medium.  This  figure  is  a  series  of  snapshots 
of  the  waveform  at  the  ranges  indicated  on  the  margin  of  the 
figure.  At  0=6.45,  and  again  at  0=14.3,  the  fundamental 
is  nearly  completely  regenerated.  Figure  18  presents  the 
amplitudes  of  the  first  three  harmonics  as  a  function  of 
range;  it  is  apparent  that  most  of  the  energy  has  been 
returned  to  the  fundamental  at  the  recurrence  distances, 


Retarded  Time,  T 


Interaction  of  Two  Solitons,  A  and  B,  in  a 
Lossless,  Nonlinear,  and  Dispersive  Medium 


(same  scale  for 
all  graphs) 


Retarded  Time , T 


Waveforms  of  a  Monofrequency  Signal 
Propagating  in  a  Lossless,  Nonlinear,  and 
Dispersive  Medium,  at  Selected  Ranges. 


Scaled  Range,  a 


Harmonic  Amplitudes  as  a  Function  of  Rang 
for  a  Monofrequency  Signal  Propagating  in 
a  Lossless,  Nonlinear,  and  Dispersive 
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by  the  coupled  effects  of  nonlinearity  and  dispersion. 

Weak-Fini t e- Am pi i tude  Propagation 
of  Waves  from  Finite  Sources 

The  numerical  method  is  capable  of  simulating  weak- 
finite-amplitude  propagation  from  finite  sources,  by  an 
appropriate  choice  of  parameters.  If  the  Rayleigh  distance 
(or  plane-wave  collimation  distance)  is  given  by 

r  =  A/ X  ,  (4.2) 

o  o 

then  the  scaled  Rayleigh  distance  is 

a  =  8ekr  .  (4.3) 

o  o 

This  quantity  is  also  called  the  scaled  source  level;  as  the 
quantity  kr^  is  an  indication  of  the  prominence  of  spherical 
spreading,  a  large  value  of  Oq  implies  that  the  nonlinear 
interaction  will  be  qua s i- plana r ,  and  a  small  value  of  Oq 
implies  that  the  nonlinear  interaction  will  extend  into 
the  spherical  spreading  region,  at  least  if  the  signal  level 
is  not  reduced  too  far  by  spreading  or  viscous  losses. 
Weak-finite-amplitude  propagation  thus  may  be  modeled  by 
using  small  values  of  cr  . 

The  scaled  Rayleigh  distance  is  used  to  define  a  new 
scaled  range  R .  If  r  is  the  physical  range  and  a  is  the 
scaled  range  as  previously  defined,  then  R  =  r/rQ  =  o/n  . 

The  acoustic  Reynold's  number  T  and  the  attenuation 


v  ^  f'  I  SMxXi 


r  "i '  V" 


so  that 

r=Oo/aQ.  (4.6) 

The  parameter  aQ  defines  the  amount  of  absorption  which 
occurs  within  the  Rayleigh  distance. 

The  combined-primary-wave  attenuation  coefficient  is 
def ined 

aT  =  ai  +  a2  -  a_  (4.7) 


so  that  the  absorption  within  the  Rayleigh  distance  is 


aT  =  ai  +  az  -  a_  . 


(4.8) 


The  results  published  by  Fenlon  and  McKendree  (1979) 
are  normalized  with  respect  to  the  parametric  array 
response  given  by  Westervelt  (1963),  which  may  be  written 


P_W  (R)  =  (P*/2Y^aTR)e~a-R, 


(4.9) 


where  R  is  the  mea n-ca r r i e r - f r e q uency  range,  P*  is  a 
reference  pressure,  YQ  is  the  frequency  downshift  ratio, 


is  the  combined-primary-wave  attenuation  coefficient,  and  a_ 

is  the  difference-frequency  attenuation  coefficient  within 

the  carrier's  Rayleigh  distance.  Thus,  to  convert  the 

parametric  gain  function  values  of  Fenlon  and  McKendree  to 

pressure,  it  is  necessary  to  multiply  the  parametric  gain 

functions  by  2a„YzR  ,  bearing  in  mind  that  the  numerical 
i  o 

solution  uses  R  to  denote  the  reference-frequency  range. 

An  additional  factor  of  4  is  introduced  by  a  difference  in 
the  definition  of  the  reference  pressure 

P*  =3PiP2koro/PoCo’  <4-10> 

as  the  numerical  solution  contains  the  implicit  assumption 

1 

that  the  peak  scaled  pressure  is  unity,  so  that  PiP2  =  4 
for  a  bifrequency  signal.  The  factor  lorQ  is  inversely 
proportional  to  the  square  of  the  frequency;  thus,  the 
actual  pressure  level  predicted  by  the  weak-finite-amplitude 
gain  function  is 

p!FA(R)  =  Gq(R)  -20  log10p(R)  -20  log  j  „  yJmJ(  4.  M ) 

Other  beam  shapes  than  Bessel  beams  from  a  plane 
piston  projector  require  different  definitions  of  the 
Rayleigh  distance,  and  corresponding  differences  in  the 
parameters  of  the  numerical  solution;  for  example,  in  the 
case  of  a  Gaussian  beam,  the  ranges  are  one-half  those  for  a 
plane  piston  projector,  and  the  required  attenuation 


coefficients  are  doubled. 


Figures  19,  20,  and  21  show  the  difference-frequency 
levels  as  a  function  of  range,  as  predicted  by  the  weak- 
finite-amplitude  solution  in  solid  curves,  and  as  predicted 
by  the  numerical  solution  as  discrete  points.  Each  graph  is 
for  a  particular  frequency  downshift  ratio  and  shows  two 
curves  for  two  different  amounts  of  absorption.  Figure  22 
shows  the  difference-frequency  level  on  axis  as  a  function 
of  range  for  a  frequency  downshift  ratio  of  5.5  and  a  norma¬ 
lized  primary-wave  attenuation  coefficient  of  1.0,  which 
represents  very  high  attenuation.  These  figures  indicate 
that  large  amounts  of  thermoviscous  absorption  reduce  the 
difference-frequency  level  in  the  far  field.  The  numerical 
solution  predicts  somewhat  higher  values  in  the  far  field 
than  does  the  weak-finite-amplitude  solution;  however,  the 
latter  has  been  found  to  underestimate  the  difference-fre¬ 
quency  level  slightly  when  it  is  used  to  simulate  experi¬ 
mental  results,  as  shown  by  Fenlon  and  HcKendree  (1979). 

Simulations  of  Experiments  at 
Finite  Amplitudes 

The  experiment  performed  by  Shooter,  Muir,  and 
Blackstock  (1974)  used  a  projector  0.0762  m  in  diameter 
operating  at  a  frequency  of  454  kHz  in  fresh  water.  The 
maximum  peak  source  level  attained  was  235  dB  re  1  MPa  at 
1  meter,  corresponding  to  1  kiloWatt  peak  input  power. 

The  following  paragraph  gives  an  example  of  the  procedure 
for  conversion  from  physical  to  scaled  parameters. 
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Figure  19.  Difference-Frequency  Level  On  Axis  as 
Function  of  Range  for  a  Bifrequency  S 
Having  a  Frequency  Downshift  Ratio  of 
for  Selected  Values  of  . 
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Figure  20.  Difference-Frequency  Level  On  Axis  as 
Function  of  Range  for  a  Bifrequency  Si 
Having  a  Frequency  Downshift  Ratio  of 
for  Selected  Values  of  a^. 
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Figure  21.  Difference-Frequency  Level  On  Ax 
Function  of  Range  for  a  Bifreque 
Having  a  Frequency  Downshift  Rat 
for  Selected  Values  of  a^. 


For  a  4  54- kHz  signal  radiating  from  a  projector 


0.07b2  m  in  diameter,  the  wave  number  is 

k  =  oj  /c  =  2. 852*  106/  1480  =  1927.0  ,  (4  .  12) 

0  0  0 

and  the  Rayleigh  distance  is 

r  =  A/A  =  1.398m  .  (4.13) 

o  o 

The  source  level  is  defined 


SL  =  20  log  (r  p  ), 
O  1  0  o  o 


(4.14) 


hence,  the  peak  signal  pressure  is 


-1  SL  /20 

p  =  r  10 
r  o  o 


(4.15) 


so  that,  for  a  peak  source  level  of  235  dB  re  1  UPa  at  1 
meter  the  peak  acoustic  pressure  has  the  value 


1.  398  1  102  3  5  /2  °  =  4.02  •  1011  liPa  =  4.02  •  106  yBar,(4.l6) 


and  the  Mach  number  has  the  value 

e  =  p  /p  c  2  =  4.02  •  10  6  / 2  .  2  5  •  1014  =  1.84  •  10~4.  (4.17) 
o  o  o  o 

For  fresh  water,  the  scaled  source  level  corresponding  to 
SLq  =235  dB  in  this  experiment  thus  is 

’  =  3r  k  r  =  1.74  .  (4.18) 

o  O  o  o 

;.'m  scaled  attenuation  parameter, a^  ,  is  the  amount  of  a  fa¬ 
it  ion  which  is  suffered  by  the  reference  frequency  in 
.  ,  i  n ^  over  the  Rayleigh  distance.  As  the  thermo- 


viscous  attenuation  at  454  kHz  in  fresh  water  is  given  by 

2  -3 

a  =  6  f  =  4. 91  •  10  ,  (4.19) 

the  scaled  attenuation  parameter  has  the  value 

a  =  ar  =  6.  86  •  10_3  (4.20) 

o  o 

The  experiment  performed  by  Shooter,  Muir,  and  Black- 
stock  has  been  simulated  by  means  of  the  numerical  solution 
with  Gaussian  beam  input.  For  this  reason,  the  results 
shown  in  the  following  figures  indicate  scaled  source  levels 
twice  those  for  a  plane  piston  projector.  Figure  23  shows 
the  level  of  the  fundamental  on  axis  as  a  function  of  range, 
in  solid  curves  for  the  numerical  results  and  as  discrete 
points  for  those  ranges  at  which  experimental  data  were 
taken.  The  numerical  solution  is  slightly  low  at  small 
ranges  and  slightly  high  at  large  ranges,  as  a  Gaussian 
beam  was  used  to  approximate  the  uniform  source  excitation. 
The  incipient  saturation  of  the  medium  is  evident,  as  the 
10  dB  increase  in  source  level  from  215  to  225  dB  re  1  p  Pa 
has  a  considerably  larger  effect  on  the  far-field  sound 
pressure  level  than  the  similar  increase  from  225  to  235  dB. 
Further  increases  in  source  level  would  raise  the  level  in 
the  near  field,  i.e.  R<1,  but  would  have  a  minimal  effect  in 
the  far  field,  R >10. 

Figure  24  presents  the  beam  widths  of  the  fundamental 
to  the  -10  dB  point  as  a  function  of  range,  for  several 
scaled  source  levels.  Finite-amplitude  losses  at  the  center 


Beam  Width 


Scaled  RangeRe 1  a t i ve  to  the  Rayleigh  Distance,  R 


Figure  24.  Numerical  Simulation  of  the  Experiment  of 
Shooter,  Muir,  and  Blackstock:  Beam  Width 
of  the  454-kllz  Fundamental  to  the  -10  dB 
Point  as  a  Function  of  Range  for  Selected 
Source  Levels. 
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of  the  beam  lower  its  level  relative  to  the  edges,  thus  pro¬ 
gressively  broadening  the  beam  as  the  scaled  source  level 
increase  s . 

Figure  25  shows  the  levels  of  the  second  harmonic  as  a 
function  of  range  for  selected  scaled  source  levels,  on 
axis  and  at  a  radial  distance  equal  to  the  effective  spot 
size  of  the  Gaussian  beam,  that  is,  to  the  original  1/e 
width  at  £=0.  Figure  26  shows  the  10-dB  beam  widths  of  the 
second  harmonic  as  a  function  of  range  for  several  scaled 
source  levels. 

Bifrequency  experiments  may  be  simulated  using  the 
numerical  method  by  choosing  the  reference  frequency,  to 
which  the  equations  are  normalized,  so  that  each  carrier  and 
the  difference  frequency  lie  at  one  of  the  Fourier  eigen- 
frequencies;  that  is,  so  that  each  carrier  and  the  dif¬ 
ference  frequency  are  periodic  within  the  set  of  waveform 
samples.  The  scaled  source  level  and  other  parameters  are 
determined  in  the  manner  outlined  previously,  but  all  must  be 
specified  in  terms  of  the  reference  frequency.  Thus,  if  Ni  , 
N2,  and  K.  are  the  Fourier  eigen  frequency  numbers  which 
correspond  to  the  carriers  and  the  difference  frequency, 

then  the  nominal  attenuation  coefficient  a  as  used  by 

T 

Fenlon  and  McKendree  (  19  7  9  )  is  converted  to  t  lie  reference 

coefficient  a  by  the  expression 
o 

a  =  a  / ( N 2  +  N 2  -  N2  (A'21) 

o  T  1  2  - )  . 


The  scaled  source  level  is  influenced  by  the  frequency 
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Scaled  Range  Relative  to  the  Rayleigh  Distance,  g 


Figure  26.  Numerical  Simulation  of  the  Experiment  of 
Shooter,  Muir,  and  Blackstock:  Beam  Width 

of  the  908-kHz  Second  Harmonic  to  the  -10 
dB  Point  as  a  Function  of  Range,  for 
Selected  Source  Levels. 


dependence  of  the  Rayleigh  distance  and  the  wave  number. 
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If  a  is  the  scaled  source  level  of  the  nominal  carrier 
o 

frequency  and  peak  signal  amplitude,  then  the  reference 
scaled  source  level  is  given  by 

a*  =  ao{  j^2N_/  (Ni  +  N  2  )^J  2  }  .  (4.22) 

In  the  experiment  performed  by  Muir  and  Willette,  the 
reference  scaled  source  level  is  o*=0. 00076  and  the  norma¬ 
lized  thermoviscous  attenuation  coefficient  is  a  =0.0003. 

o 

This  experiment  used  a  projector  0.02  m  in  diameter,  oper¬ 
ating  at  418  and  482  kHz  to  generate  a  64-kHz  difference 
frequency.  For  the  numerical  solution,  N  was  chosen  as  2 
and  the  reference  frequency  as  32  kHz.  The  carriers  of  480 
and  416  kHz  corresponding  to  Nr=15  and  N2=13  are  a  good 
approximation  of  the  actual  values.  Figure  27  shows  the 
level  of  the  difference  frequency  on  axis  as  a  function  of 
range.  The  numerical  solution  is  shown  as  a  solid  curve  and 
the  values  measured  in  the  experiment  are  shown  as  discrete 
points.  Values  derived  from  the  weak-finite-amplitude 
solution  are  also  shown  for  purposes  of  comparison.  Figure 
28  presents  the  beam  widths  of  the  difference  frequency  as  a 
function  of  range  as  measured  and  as  predicted  by  the 
solutions.  The  numerical  results  are  fairly  close  to  the  ex¬ 
perimental  values.  Cases  such  as  this,  involving  large  fre¬ 
quency  downshift  ratios,  very  small  scaled  source  levels, 
and  large  ranges,  are  relatively  difficult  for  the  numerical 
solution  to  simulate.  For  such  cases,  the  weak-finite-ampli- 
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28.  Numerical  Simulation  of  the  Experiment  of 
Muir  and  Willette:  Beam  Widths  of  the 
Difference  Frequency  as  a  Function  of  Range. 


tude  solution  is  computationally  more  efficient. 

The  experiment  performed  by  Eller  (1974)  used  a  pro¬ 
jector  0.02  m  in  diameter,  operating  at  a  mean  carrier  fre¬ 
quency  of  1435  kHz  in  fresh  water.  The  numerical  simulation 
used  N  =29,  =28,  and  N_  =1,  with  a  reference  frequency  of 

50  kHz,  equal  in  this  case  to  the  difference  frequency. 
Figure  29  shows  the  difference-frequency  level  on  axis  as  a 
function  of  range,  as  measured  by  Eller  and  as  predicted  by 
the  numerical  and  weak-finite-amplitude  solutions.  Figure 
30  shows  the  beam  widths  to  the  -10  dB  point  as  a  function 
of  range,  as  predicted  and  as  measured. 

Numerical  Simulations  at  Strong 
Finite  Amplitudes 

The  fundamental  assumption  of  the  weak-finite-amplitude 

solution  is  that  the  carriers  are  not  subject  to  finite- 

amplitude  losses.  As  the  scaled  source  level  is  increased, 

finite-amplitude  losses  in  the  carriers  become  significant. 

The  reduced  carrier  levels  in  turn  reduce  the  amount  of 

energy  which  may  be  downshifted  to  the  difference  frequency. 

As  the  scaled  source  level  becomes  very  large,  the  nonlinear 

interaction  is  restricted  to  the  plane-wave  collimat  ion  zone 

within  the  Rayleigh  distance. 

The  fundamental  parameters  of  the  scaled  nonlinear 

propagation  problem  are  the  scaled  source  level  o  ,  the 

o 

attenuation  within  the  Rayleigh  distance  of  the  reference 
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Figure  30.  Numerical  Simulation  of  the  Experiment  of 
Eller:  Beam  Widths  of  the  Difference 

Frequency  as  a  Function  of  Range. 
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frequency  aQ  ,  and  the  nature  of  the  signal  and  projector. 

A  bifrequency  pair  may  be  characterized  by  its  frequency 

downshift  ratio  Y  . 

o 

The  following  numerical  simulations  are  of  two  types: 
numerical  simulations  of  previously  published  research,  and 
simulations  of  experiments  with  scaled  parameters  chosen  so 
as  to  emphasize  the  effect  of  one,  or  to  combine  the  effects 
of  several. 

The  scaled  source  level  indicates  the  strength  of  the 
nonlinear  interaction  relative  to  other  effects  of  propa¬ 
gation.  The  choice  of  scaled  source  level  values  at 

l  /  2 

approximately  half-decade  increments,  or  powers  of  10  ,  as 

10,  3,  1,  0.3,  0.1,  and  so  on,  permits  a  reasonable  view  of 
the  continuum  of  scaled  source  levels.  The  attenuation 
coefficients  aT  are  chosen  in  decimal  increments  to  repre¬ 
sent  high  attenuation,  which  dominates  the  nonlinearity,  and 
lower  attenuation,  which  is  dominated  by  nonlinearity. 

As  the  Fourier  harmonics  are  have  integer  ratios,  integer 
pairs  are  chosen  to  represent  a  bifrequency  signal,  e.  g. 
N^=6  and  =5  for  a  frequency  downshift  ratio  of  3.3.  Thus, 
frequency  downshift  ratios  of  5.5  and  10.5  are  used  in  the 
following  simulations. 

Figure  31  shows  the  difference-frequency  level  on  axis 
as  a  function  of  range,  as  predicted  by  the  numerical  solu¬ 
tion  for  various  scaled  source  levels  and  as  given  by  the 
weak-finite-amplitude  solution  for  a  frequency  downshift 
ratio  of  5.5.  Figure  32  shows  the  same  functions  for  a  fre- 
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quency  downshift  ratio  of  10.5.  For  each  value  of  the 
scaled  source  level,  the  difference-frequency  level  is  the 
same  as  the  weak-finite-amplitude  solution  at  small  ranges. 
As  the  range  increases,  saturation  of  the  carriers  causes 
the  difference  frequency  level  to  decrease.  The  greater 
the  scaled  source  level,  the  sooner  the  departure  occurs. 

For  a  scaled  source  level  of  1  or  greater,  the  increased 
source  level  is  almost  entirely  canceled  by  finite-amplitude 
losses;  the  nonlinear  interaction  is  saturated  within  the 
Rayleigh  distance. 

Figure  33  shows  the  difference-frequency  level  on  axis 
as  a  function  of  range,  as  predicted  by  the  numerical  sol¬ 
ution  and  as  given  by  the  weak-finite-amplitude  solution,  at 
various  scaled  source  levels  for  a  frequency  downshift  ratio 
of  5.5  and  a  reference  combined-primary-wave  attenuation 
coefficient  of  1.0.  In  this  instance  t he rraoviscous  losses 
are  large  enough  to  reduce  the  carrier  levels  to  a  large 
extent  within  the  critical  range. 

The  on-axis  pressure  in  a  Gaussian  beam  propagating  in 
a  linear  manner  obeys  the  following  function  of  range: 

|PU(»)I“  |jpa)(0)i/(l+R2)  /ZJe  (4.23) 

If  the  amplitude  of  the  relative  difference-frequency  level 
derived  from  numerical  results  is  divided  by  the  right-hand 
side  of  Equation  (4.23)  the  effects  of  thermoviscous  and 
spreading  losses  are  canceled  in  that  region  of  the  far 
field  wherein  only  the  linear  loss  mechanisms  represented 
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in  Equation  (4.23)  remain.  Figures  34,  35,  and  36  each 
show  the  difference-frequency  level  adjusted  for  thermo- 
viscous  and  spreading  losses  as  a  function  of  range,  for 
various  frequency  downshift  ratios  and  selected  values  of 
reference  attenuation  coefficient.  These  may  be  called 
saturation  curves,  as  the  saturation  of  the  nonlinear 
interactions  is  reflected  in  their  each  attaining  a  constant 
value.  From  such  curves,  the  extrapolated  source  level  may 
be  deduced  by  finding  the  saturation  level  and  referring  it 
back  to  the  source. 

The  modified  Bessel-Fubini  solution  derived  by  Fenlon 
(1972)  provides  a  good  model  of  bifrequency  interaction  in 
lossless  media  near  the  source.  Figures  37  and  38  show 
predicted  pressures  in  a  bifrequency  plane  wave  due  to  the 
source  terms  of  the  Bessel-Fubini  solution,  and  as  derived 
from  the  numerical  solution.  The  spherically-spreading 

results  for  a  scaled  source  level  o  =20  are  also  shown. 

o 

These  are  very  similar  to  the  plane-wave  results. 

The  numerical  work  of  Bakhvalov  et  al.  (1978,1979), 
concerning  initially  monotonic  waves  in  t he rmo vis cous  media 
at  strong  finite  amplitudes,  has  been  simulated  numerically. 
Figure  39  shows  the  waveforms  at  R=0.09  and  several  values 
of  the  radial  coordinate,  for  a  scaled  source  level  of  20 
and  a  reference  attenuation  coefficient  of  0.03.  The  cor¬ 
responding  figure  from  Bakhvalov's  results  is  included  for 
purposes  of  comparison.  Figure  40  shows  the  waveform  on 
axis  at  several  ranges,  for  a  scaled  source  level  of  10, 
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Figure  35.  Difference-Frequency  Level  On  Axi 
Compensated  for  Spherical  Spreadi 
The rmo v i s c ou s  Losses,  With  a  Free 
Downshift  Ratio  of  10.5  and  aT=4. 
for  Selected  Scaled  Source  Levels 
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an  attenuation  coefficient  of  0.03,  and  a  beam  shape  at  the 
source  defined  by  1-  £ 4  to  £=1  and  zero  beyond  £=1.  Here 
also  a  comparison  figure  is  included. 

Bakhvalov  et  al.  included  in  their  results  figures 
illustrating  the  region  of  shock  formation,  specified  as  a 
"quasi-discont inuous  wave",  as  a  function  of  range  and 
radial  distance  for  various  scaled  source  levels  and  amounts 
of  absorption.  Figure  41  shows  regions  of  shock  formation, 
defined  as  that  region  in  which  the  harmonics  equal  or  exceed 
threshold  values  given  on  the  figure.  The  corresponding 
figure  from  Bakhvalov's  results  is  included  for  comparison. 

The  numerical  results  are  based  upon  the  attainment  of 
specified  harmonic  levels,  rather  than  on  the  slope  of  the 
discontinuity  as  was  used  by  Bakhvalov  et  al.;  therefore  the 
appearance  of  the  former  and  latter  is  not  precisely  the  same. 
A  waveform  having  50  percent  second  harmonic,  33  percent  third 
harmonic,  and  so  on  in  the  correct  relative  phase,  will  have 
a  discontinuity  extending  from  the  positive  peak  to  the 
negative  peak  of  the  waveform. 
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Chapter  5 

CONCLUSIONS 

The  first  section  of  the  conclusions  deals  with  the 
results  of  tests  of  the  numerical  solution  to  determine  the 
accuracy  which  is  exhibited  in  these  results.  In  the  next 
section,  the  limits  of  applicability  of  the  weak-finite- 
amplitude  theory  are  discussed.  The  topic  of  near-field 
calibration  of  parametric  sources  is  then  considered.  The 
conversion  efficiencies  which  may  be  expected  of  parametric 
arrays  is  indicated  by  analysis  of  the  results  of 
numerical  simulations.  Finally,  suggestions  are  given  for 
further  research  in  the  numerical  solutions  of  nonlinear 
wave  equations. 

Tests  of  the  Numerical  Solution 

Table  5.1  shows  the  comparison  between  the  matched 
asymptotic  lossless  initially  monotonic  plane-wave  solution 
due  to  Blackstock  and  the  results  of  numerical  solution  of 
the  lossless  Burger's  equation.  At  each  of  several  ranges 
and  frequencies,  the  difference  between  the  numerical 
solution  and  Blackstock's  solution  is  presented  as  a 
percentage  of  the  latter.  The  harmonic  levels  of  the 
numerical  solution  are  within  a  few  percent  of  Blackstock's 


solution  throughout  the  range.  This  degree  of  accuracy 
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Table  5.1.  Numerical  Harmonic  Levels  of  the  First 

Three  Harmonics  of  a  Monofrequency  Signal  at 
Selected  Ranges,  and  the  Values  of  the  Matched 
Asymptotic  Solution  due  to  Blackstock  (1966). 


SIGMA 

NUMERICAL 

ANALYTIC 

(a) 

B  (a  ) 

1 

ERROR  IN 
PERCENT 


0.5 

.9649 

.  969  1 

-0.4 

1.0 

.  8747 

.8801 

-0.6 

2.0 

.  6435 

.6484 

-0.8 

3.0 

.4891 

.4937 

-0 . 9 

5.0 

.  3274 

.  3333 

-1  .  7 

(b) 

B  (a  ) 

2 

0.5 

.  2274 

.  2298 

-1 . 0 

1 .0 

.  3483 

.3528 

-1  .  3 

2.0 

.3102 

.3134 

-1 . 0 

3.0 

.  2444 

.2428 

0.7 

5.0 

.  1682 

.1667 

0.9 

(c) 

B  (a  ) 

3 

0.5 

.0791 

.0813 

-2.7 

1  .0 

.2102 

.  2060 

2.0 

2.0 

.  2025 

.2070 

-2 . 2 

3.0 

.  1581 

.1610 

-1  .  8 

5.0 

.  1055 

.1111 

-5.0 

- « 


confirms  the  validity  of  the  numerical  solution  of  the 
lossless  Burger's  equation. 

Tables  5.2a,  b,  and  c  show  the  far-field  difference- 
frequency  level  times  expfo  / T) ,  as  predicted  by  the  solution 
of  Fenlon  (1972)  and  as  derived  from  the  numerical  solution 
of  the  lossy  Burger's  equation.  The  agreement  between 
analytic  and  numerical  results  is  fairly  good.  For  some  of 
these  simulations,  as  many  as  4000  iterations  of  the 
nonlinear  operator  were  required.  The  presence  of  bias  in 
the  results  is  therefore  to  be  expected,  though  the  largest 
error  does  not  exceed  12  percent. 

The  numerical  simulation  of  bifrequency  propagation  in 
a  lossy  medium  requires  the  numerical  solution  to  apply  the 
nonlinear  operator  simultaneously  at  several  frequencies,  so 
as  to  simulate  s e 1 f -con volu t ion  of  the  input  signal,  and  to 
apply  the  correct  thermoviscous  losses  to  each  of  the 
resulting  frequency  components.  In  fact,  the  analytic 
far-field  difference-frequency  levels  are  reproduced  by  the 
numerical  solution  with  an  amount  of  error  not  much  larger 
than  may  be  expected  due  to  the  bias  of  the  lossless 
nonlinear  operator  itself.  It  is  evident  that  the  number 
of  waveform  samples  and  the  number  of  steps  per  unit  range 
will  not  significantly  affect  the  results,  so  long  as  each 
is  large  enough.  The  ability  of  the  numerical  solution  to 
simulate  bifrequency  propagation  in  a  lossy  medium 
indicates  that  both  nonlinearity  and  viscous  losses  are 


properly  modeled. 
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Table  5.2.  Difference-Frequency  Level  in  the  Far  Field 
Times  exp(a/r)  As  Predicted  by  the  Numerical 
Solution  and  by  the  Solution  due  to  Fenlon 
(1972),  for  Selected  Frequency  Downshift  Ratios 
and  Values  of  T. 


GAMMA 

ANALYTIC 
(Fenlon, 1972) 

NUMERICAL 

ERROR  IN 
PERCENT 

(a)  frequency 

downshift 

ratio  5.5 

11 

.0428 

.043 

-0.5 

16 

.0583 

.057 

-2 . 2 

22 

.0724 

.070 

-3.3 

31 

.0847 

.081 

-4.3 

44 

.0885 

.085 

-3.9 

62 

.0817 

.077 

-5.7 

88 

.0678 

.065 

-4  .  1 

124 

.0530 

.051 

-3.8 

(b )  frequency 

downs h i f  t 

ratio  7.5 

15 

.0314 

.031 

-1  .  3 

21 

.0415 

.040 

-3.6 

30 

.0531 

.051 

-4.0 

42 

.0620 

.060 

-3.2 

60 

.0649 

.060 

-7.6 

85 

.0598 

.056 

-6.4 

120 

.0497 

.047 

-5.4 

170 

.0387 

.038 

-1.8 

(c)  frequency 

downsh i f  t 

ratio  10.5 

21 

.0224 

.022 

-1.8 

30 

.0301 

.029 

-3.7 

42 

.0380 

.036 

-5.3 

59 

.0443 

.040 

-9.7 

84 

.0464 

.041 

-11.6 

119 

.0427 

.041 

-11.0 

168 

.0355 

.032 

-9.8 

238 

.0276 

.025 

-9.4 

)• 

I 

f 
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Simulated  propagation  in  a  lossless,  dispersive  medium 
is  a  severe  test  of  the  numerical  solution.  Since 
reasonably  accurate  simulations  of  dispersive  propagation 
result  from  the  use  of  the  numerical  solution,  it  may  be 
assumed  that  the  method  of  solution  is  valid. 

The  preceding  tests  apply  to  plane-wave  propagation, 
for  which  the  pertinent  equation  is  equation  4.23.  In  order 
to  accommodate  spreading  as  in  equation  4.22,  the  diffrac¬ 
tion  term  must  be  handled  by  the  method  of  diffusion-in¬ 
frequency.  The  present  numerical  implementation  works  well 
out  to  R  on  the  order  of  20;  spherical  spreading  may  be 
assumed  to  prevail  beyond  this  range.  The  numerical 
solution  is  designed  to  allow  the  user  to  choose  a  value 
of  R  beyond  which  spherical  spreading  is  assumed.  If  this 
value  is  chosen  carefully,  the  resulting  error  is  not 
intolerable . 

The  ability  of  the  numerical  solution  to  produce 
correct  results  when  a  combination  of  factors  is  present, 
such  as  nonlinearity  and  attenuation  or  nonlinearity  and 
dispersion,  indicates  that  the  implementation  which  was 
used,  regarding  each  of  these  effects  as  independent  over 
small  distances,  is  sufficiently  accurate.  There  exist 
certain  areas  in  which  the  computational  efficiency  could 
be  improved.  These  are  discussed  in  the  suggestions  for 
further  research.  The  remaining  computational  problems  are 
also  discussed  in  the  last  section. 


Limits  of  Applicability  of  the 
Weak-Fini te-Amplitude  Theory 

The  results  of  numerical  analysis  confirm  that  the  weak- 
finite-amplitude  solution  is  invalidated  by  f i ni te-amp 1 it ude 
losses  in  the  carriers.  For  relatively  lossless  cases,  a 
simple  constraint  suffices  to  ensure  the  applicability  of 
the  weak-finite-amplitude  theory:  the  maximum  range  of 
interest  must  be  considerably  less  than  the  critical  range 
of  either  of  the  carriers,  so  that  they  do  not  lose  a 
significant  amount  of  energy  in  generating  harmonics  or 
modulation  frequencies. 

The  limiting  case  as  the  scaled  source  level  becomes 
infinite,  physically  the  case  for  very  high  frequencies 
or  very  large  apertures,  is  that  of  collimated  plane-wave 
propagation;  the  nonlinear  interaction  is  confined  to  the 
very  near  field  of  the  source.  The  largest  scaled  source 
levels  actually  used  in  computation,  on  the  order  of  20, 
are  virtually  indistinguishable  from  plane-wave  propagation 
in  the  near  field.  In  these  cases,  the  carriers  of  a 
bifrequency  signal  do  not  go  into  shock,  but  attain  second 
harmonic  levels  of  20  to  25  percent  at  most.  The  greatest 
amount  of  energy  is  transferred  to  the  sum  frequency  and  its 
harmonics,  which  will  attain  a  shock  profile  in  a  relatively 
lossless  medium. 

In  very  lossy  media,  the  wea k-f i n i t e -amp 1 i t u d e  solution 
is  applicable  at  higher  scaled  source  levels  than  in 
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less  lossy  media.  The  linear  losses  reduce  the  carrier 
amplitudes  within  the  near  field,  thus  suppressing  part  of 
the  nonlinear  interaction.  For  example,  with  a  frequency 
downshift  ratio  of  5.5  and  a  scaled  source  level  of  0.1, 
the  numerical  and  weak-finite-amplitude  results  differ  by 
8  dB  at  R=20  if  the  combined-primary-wave  attenuation 
coefficient  is  0.001;  whereas  if  it  is  increased  to  1.0, 
the  numerical  and  weak-finite-amplitude  results  are 
virtually  the  same  out  to  R  =  30 . 

Near-Field  Behavior  and  Calibration 
of  Parametric  Sources 

The  numerical  solution  provides  a  means  of  modeling  the 
operation  of  parametric  arrays  in  both  the  unsaturated  and 
saturated  regimes.  The  pressure  waveform  or  spectrum  may 
be  obtained  from  the  numerical  solution  at  any  point 
within  the  Fresnel  zone,  as  the  solution  of  the  diffraction 
term  is  based  upon  a  Fresnel  approximation.  Zemanek  (1971) 
published  a  paper  which  investigates  the  near-field 
behavior  of  the  beam  from  an  axisymmetric  vibrating  piston 
of  finite  extent.  His  method  has  been  implemented  with 
excellent  results,  but  as  it  is  vastly  slower  than  the 
diffusion-in-frequency  technique,  it  is  not  suited  for 
iterative  computations. 

The  effective  length  of  a  parametric  source  operating 
in  a  relatively  lossless  medium  may  be  much  longer  than  the 
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Rayleigh  distance  of  the  projector  at  the  difference 

frequency,  which  in  many  practical  cases  may  be  too  long  to 

allow  calibration  of  the  difference-frequency  signal  in  an 

anechoic  facility  of  practical  size.  The  weak-finite- 

amplitude  theory  is  capable  of  predicting  the  near-field 

behavior  in  the  Fresnel  zone  of  unsaturated  arrays.  The 
numerical  solution  as  presently  implemented  is  also 

restricted  to  the  Fresnel  zone,  but  will  accommodate 

saturated  parametric  arrays.  Thus  calibration  of  the 

difference-frequency  characteristics  of  a  saturated 

parametric  source  is  possible  within  its  near  field  by 

comparison  with  the  results  of  the  numerical  solution. 

Conversion  Efficiencies  of 
Parametric  Arrays 

The  amplitudes  of  carriers  from  a  parametric  source  are 
reduced  by  finite-amplitude  losses,  by  spreading,  and  by 
therraoviscous  losses.  Eventually,  they  will  be  reduced  to  a 
point  where  they  no  longer  transfer  a  significant  amount  of 
energy  to  the  difference-frequency  signal.  The  difference 
frequency  will  propagate  thereafter  under  the  influence  of 
spreading  and  thermovis cous  losses.  From  the  difference- 
frequency  level  in  this  region,  an  extrapolated  source  level 
may  be  deduced.  A  d i f f e r e nce-f r equen cy  signal  radiated 
from  the  source  position  at  this  level  would  give  rise  to 
the  pressure  observed  in  the  far  field  as  a  result  of  the 
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parametric  interaction. 

It  is  of  interest  to  determine  what  may  be  the  extra¬ 
polated  source  level  of  the  difference  frequency  relative 
to  the  level  of  the  carriers.  This  may  be  done  by  analysis 
of  the  numerical  results.  Values  of  absorption,  scaled 
source  level,  frequency  downshift  ratio  or  modulation  type, 
and  beam  shape  are  chosen.  Tie  signal  is  allowed  to 
propagate  until  nonlinear  energization  of  the  difference 
frequency  is  seen  to  have  ended.  The  far-field  difference- 
frequency  pressure  is  converted  to  a  pressure  at  the  source 
by  application  of  Equation  (4.12),  or  its  appropriate 
analogue  for  other  beam  types,  and  this  equivalent  source 
pressure  is  expressed  in  dB  relative  to  the  peak,  pressure 
which  was  actually  transmitted. 

The  numerical  solution  requires  the  peak  pressure  at  the 
source  to  be  unity  to  account  for  the  correct  scaling  of 
the  nonlinear  operator.  Thus,  the  numerical  value  obtained 
from  the  program  at  any  frequency,  range,  and  radial 
distance  is  that  which  would  be  produced  by  a  unit  peak 
source.  The  Rayleigh  distance  does  not  appear  explicitly  in 
this  calculation.  The  following  figures  show  the  extra¬ 
polated  source  levels  of  the  difference  frequency  as  a 
function  of  scaled  source  level  of  the  carriers  for  selected 
combined-primary-wave  attenuation  coefficients. 

Figure  42  shows  the  difference-frequency  extrapolated 
source  levels  for  a  bifrequency  signal  having  a  frequency 
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downshift  ratio  of  5.5,  and  Figure  43  for  a  frequency 
downshift  ratio  of  10.5.  As  an  example  of  their  use,  con¬ 
sider  a  bifrequency  signal  having  a  frequency  downshift 
ratio,  y  ,  of  5.5,  a  scaled  source  level,  Oq,  of  0.3,  and 
a  combined-primary-wave  attenuation  coefficient,  aT,  of  1.0. 
By  reference  to  Figure  42  it  may  be  determined  that  the 
extrapolated  source  level  of  the  difference  frequency  is 
38  dB  below  the  peak  level  of  either  of  the  carriers  at  the 
source.  The  numerical  solution  may  be  used  as  shown  above  to 
generate  families  of  curves  of  this  kind  for  any  modulation 
type  and  range  of  source  levels,  limited  only  by  the  amount 
of  space  available  for  storage  of  waveform  samples,  and  the 
amount  of  time  available  for  computation. 

Suggestions  for  Further  Research 

The  method  of  numerical  solution  employed  in  this  thesis 
produces  good  results  in  most  cases.  The  validity  of  the 
basic  approach  to  the  problem  is  confirmed  by  the  quality 
of  its  results.  The  method  is  subject  to  failure  in  certain 
respects,  and  as  these  impose  the  major  restrictions  on  the 
usefulness  of  the  numerical  solution,  it  is  the  purpose  of 
this  section  to  discuss  the  causes  of  the  computational 
problems  and  the  benefits  which  may  be  expected  from 
alleviating  them. 

The  simple  linear  interpolator  used  in  the  solution  of 
the  lossless  Burger's  equation  is  highly  satisfactory,  but 
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it  requires  a  very  small  step  size,  Ao  ,  if  t he  highest 
input  frequency  is  a  large  multiple  of  the  reference 
frequency  or  Fourier  fundamental.  If  in  addition  the  scaled 
source  level,  0Q,  is  large,  the  number  of  iterations 
needed  to  reach  the  far  field,  R>>1,  can  become  very  large. 
The  higher-quality  interpolators  which  were  tried  in  solu¬ 
tions  to  the  lossless  Burger's  equation  provided  the  same 
accuracy  as  the  linear  interpolator  with  a  much  larger  step 
size  until  a  discontinuity  was  formed,  at  which  point 
instabilities  would  develop  in  the  vicinity  of  the  discon¬ 
tinuity.  A  method  using  a  high-quality  interpolator  in 
regions  where  a  discontinuity  will  not  cause  difficulties 
could  improve  the  efficiency  of  the  numerical  solution. 

The  diffusion-in-frequency  technique  is  implemented  by 
means  of  an  implicit  Cr ank-Nic hols  on  formula.  This  method 
is  stable  and  efficient,  and  works  well  out  to  R  on  the 
order  of  20,  at  which  point  the  error  in  the  level  on  axis 
is  0.3  dB  and  the  largest  errors  in  bearawidth  are  about  10 
percent.  The  reason  for  the  errors,  which  become  very 
large  at  larger  values  of  R,  has  been  found  to  be  the  way 
in  which  the  beam  profile  is  sampled  rather  than  the  Crank- 
Nicholson  formula  itself.  Other  methods  of  solution,  which 
are  less  efficient  than  the  Crank-Nicholson  technique  but 
used  the  same  beam  sampling  technique,  have  been  tried  and 
found  subject  to  the  same  sampling  problem.  This  problem  is 
associated  with  the  spreading  of  the  radial  coordinate  sys¬ 


tem  as  it  follows  the  beam  propagating  away  from  the  source. 
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Spreading  is  characteristic  of  propagation  from  a  source 
of  finite  extent.  In  the  far  field  of  a  linear  projector, 
the  radial  beamwidth  is  proportional  to  the  range  in  the 
direction  of  propagation,  and  the  angular  beamwidth  is 
constant.  Therefore,  if  a  given  number  of  beam  samples  at 
a  given  spacing  suffices  near  the  source,  then  either  the 
number  of  samples  or  the  radial  distance  increment  must  be 
increased  at  larger  distances  if  the  edges  of  the  beam  are 
to  remain  visible  as  it  propagates.  The  latter  option,  to 
increase  the  radial  sample  spacing  when  necessary,  has  been 
chosen  for  two  reasons.  First,  the  amount  of  storage 
available  is  fixed,  and  thus  the  number  of  samples  cannot  be 
allowed  to  grow  without  limit,  and  second,  the  time  required 
for  calculation  increases  rapidly  as  the  number  of  radial 
samples  increases. 

The  expansion  of  the  radial  sample  spacing  is  done  in 
the  following  way.  The  beam  samples  at  the  larger  spacing 
are  interpolated  or  selected  from  among  those  presently 
existing,  which  represent  samples  at  a  smaller  sample 
spacing.  The  samples  at  expanded  coordinates  outside  the 
domain  previously  defined  are  all  assumed  to  be  zero. 

This  technique  leads  to  progressive  deterioration  of  the 
results  beyond  R=20,  with  considerable  distortion  of  the 
beam  shape  by  R=50. 

Spherical  spreading  may  be  assumed  to  apply  at  large 
ranges.  The  manner  in  which  this  option  is  implemented  in 
the  numerical  solution,  although  it  works  very  well  in 
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modeling  linear  diffraction,  leads  to  difficulties  in  the 
nonlinear  cases,  as  the  phase  distribution  of  the 
nonlinearly  generated  frequencies  cannot  be  predicted  a 
priori.  In  this  thesis,  the  Crank-Nichols on  formula  is 
used  exclusively  for  results  confined  to  R  less  than  about 
30,  and  for  longer  ranges  a  transition  zone,  e.  g.,  from 
R=10  to  R=20 ,  is  defined  and  the  Crank-Nichols on  and 
imposed-spherical-spreading  results  are  matched  at  the 
limits  of  the  transition  zone. 

Advances  in  the  nonlinear  operator  and  in  the  expan- 
ding-grid  diffraction  technique  should  allow  rapid  and 
convenient  numerical  solution  of  many  problems  which  are 
described  by  the  second-order-nonlinear  wave  equation, 
more  quickly  and  more  accurately  than  the  present 
implementation.  The  techniques  of  Zemanek  (1971)  and  of 
Sziklas  and  Siegman  (1974)  should  be  useful  points  of 
departure  in  improving  the  numerical  solution  of  the  dif¬ 
fraction  term.  Zemanek  demonstrated  the  usefulness  of  direct 
numerical  integration  of  the  diffraction  integral.  Although 
this  method  is  not  fast  enough  to  be  used  in  an  iterative 
numerical  solution,  it  might  be  used  in  performing  the  radial 
coordinate  expansion  which  is  required  from  time  to  time, 
without  the  errors  associated  with  the  present  sampling  tech¬ 
nique.  Szikias  and  Siegman  show  a  method  for  evaluating  the 
diffraction  integral  in  spherical  coordinates  by  the  use  of 
weighted  Fourier  transforms,  in  a  modified  system  of 
coordinates.  These  coordinates  are  such  that  a  spherically 
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APPENDIX 

LARGE-FINITE-AMPLITUDE  COMPUTER  PROGRAM 
SOURCE  LISTINGS 

The  following  source  listings  give  the  code  for  the 
two  programs  used  in  deriving  the  numerical  results  pre¬ 
sented  in  this  thesis.  Other  versions  differ  only  in  slight 
details  of  the  coding,  usually  in  the  format  of  output. 

The  program  named  PRPGT2  is  the  main  program  which  solves 
the  plane-wave  nonlinear  propagation  problem.  The  program 
named  PRPGT8  solves  the  case  of  axisymmetr ic  three-dimen¬ 
sional  propagation.  Both  of  these  programs  use  the  solution 
to  Burger's  equation  expressed  in  Equation  (3.4),  as  this 
was  the  formulation  used  to  derive  most  of  the  results 
presented  in  this  thesis. 


nnoonnnnnnnnnn  nnnnnnnnnnnnnnnnnnnnnonnnnnnnrinonn 
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PRPGT8  RELEASE  13.0  —  8  MARCH  1981  —  F  S  MCKENDREE 
FOR  FURTHER  INFORMATION  CONTACT - 

FRANCIS  S  MCKENDREE  ! 
APPLIED  RESEARCH  LABORATORY  . 
ROOM  1 23A  L 
P  0  BOX  30  i 
STATE  COLLEGE,  PA  16801 
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THIS  PROGRAM  SOLVES  THE  SECOND-ORDER-NONLINEAR  PARABOLIC  WAVE 
EQUATION. 


REVISION  4  FEBRUARY  1981  —  WAVEFORM  PRINTOUT  CAPABILITY  ADDED 
AT  RANGES  CALLED  FOR  IN  'SMAX'  AND  RADIAL  COORDINATES  GIVEN  IN 
ARRAY  ' EP RNT ‘ .  MAXIMUM  8  RADIAL  VALUES  AT  EACH  DISPLAY  CYCLE; 
SPECIFY  RADIAL  VALUE  LESS  THAN  ZERO  TO  END  DISPLAY  AT  GIVEN 
RANGE,  OR  SMAX  LESS  THAN  ZERO  TO  END  ITERATIONS. 

REVISION  8  MARCH  1981  --  ( 1 )  DELTA  SIGMA  MAY  BE  CHANGED  AT  THE 
BEHEST  OF  THE  USER  TO  THE  NEW  VALUE  'DSIGN'  AT  THE  RANGE 
'SIGCH' ;  (2)  AN  AXIAL  BEAM  LEVEL  IS  DEFINED  AS  THE  MAXIMUM 
FOUND  WITHIN  THE  1/2  DB  WIDTH  OF  A  LINEAR  GAUSSIAN  BEAM  AT  EACH 
RANGE,  AND  IS  DISPLAYED  AND  USED  AS  THE  ZERO-DB  LEVEL  FOR  THE 
BEAM  WIDTH  CALCULATIONS. 


COMPLEX  UO (256) ,  UZ(256),  W(256),  CF.  CFC,  SPMD(129),  AD,  AF,  SPI 
REAL  ROOTV(3)/0. 7071068.0. 5,0. 3162278/.  BW(3,129),  EPSV(4l) 

REAL  U(258,40),  ATTF(129).  £>MAG(256),  UR(2563 
REAL  QLL/1 .000/ .  RSPH/10.0/,  RDIFF/0.300/ 

INTEGER*4  TITLE!  16 )/  'PROG' ,  kAM  ','PRPC'/T8  RVELEA'/SE 

1  '13.0','  8  VMARCVH  V  19&','l  7 

INTEGER  IDAT(2).  JDAT(2) 

INTEGER  IDIFFR( 16 ) ,  ISQ(256) 

LOGICAL  NOLOSS.  NORECD,  NODIFF,  NODISP,  NOPRNT,  NOSPM 
LOGICAL  SW.  SI/ ,  Sll.  S10.  S09 ,  NONONL,  SW2 
REAL  EPRNT(8),  AXISLV(129) 


SENSE  SWITCH  OPTIONS: 

SWITCH  11  ON  TO  LIST  INITIAL  VARIABLES. 
SWITCH  10  ON  TO  SUPPRESS  DIFFRACTION. 

SWITCH  9  ON  TO  ENABLE  RECORDING  OF  DATA. 
SWITCH  8  ON  TO  SUPPRESS  LOSSES. 

SWITCH  7  ON  TO  SUPPRESS  DISPERSION. 

SWITCH  6  ON  TO  SUPPRESS  NONLINEAR  OPERATOR. 
SWITCH  5  ON  TO  READ  QLL ,  RSPH,  AND  RDIFF 
SWITCH  4  ON  TO  PRINT  INITIAL  WAVEFORMS. 


SETUP . 

FILE  'ICARD'  IS  USED  FOR  ALL  SUBSEQUENT  CONTROL  INPUT. 


READ(5,8000)  ICARD 
CALL  SSWTCH( 5 ,SW2 ) 

IF  (SW2)  READ(5,3020)  QLL,  RSPH,  RDIFF 
100  SW=. FALSE. 

200  CALL  SSWTCK(7, NODISP) 

CALL  SSWTCH(8, NOLOSS) 

NOSPM=. FALSE. 

IF  (NOLOSS  .AND.  NODISP)  NOSPM=.TRUE. 
CALL  SSWTCH( 10, NODIFF) 

CALL  SSWTCH( 6, NONONL) 

CALL  DATE(IDAT) 
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C 

C 

C 


150 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


300 

400 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


CALL  LINETT  (  TITLE  ) 

ZERO  ATTF ,  DSP,  AND  SPMD  ARRAYS 

CF=CMPLX( .0,0.) 

DO  150  1=1,129 
ATTF(I)=0. 

SPMD(I)=CF 

CONTINUE 

INPUT  CONTROL  VARIABLES: 

INL  -  NUMBER  OF  WAVEFORM/RADIAN  FREQUENCY  SAMPLES;  NUMBER  OF 
POINTS  IN  THE  NON  LINEAR  OPERATOR. 

IFB  -  NUMBER  OF  RADIAL  SAMPLE  POINTS 
ICASE  -  CASE  NUMBER  FOR  REFERENCE  PURPOSES 

NCOMP  -  NUMBER  OF  FOURIER  COMPONENTS  WHICH  ARE  TO  BE  ASSUMED  TO 
EXIST  IN  THE  INPUT  SIGNAL,  FOR  PURPOSES  OF  NORMALIZATION;  IF 
ZERO,  NO  NORMALIZATION  IS  PERFORMED. 

READ( ICARD, 8000,END=7 200)  INL,  IFB,  ICASE,  NCOMP 

IDIFFR  -  ARRAY  STORING  THE  DFT  FREQUENCIES  AT  WHICH  THE 
NEED  FOR  RADIAL  EXPANSION  IS  TO  BE  TESTED 
DETERMINE  THE  NUMBER  OF  FREQUENCIES  SELECTED  IN  'ISF'. 

READ (ICARD, 8000)  IDIFFR 
K=1 

DO  300  1=1.16 

IF  (IDIFFR(I)  .LT.  0)  GO  TO  400 
K=I 

CONTINUE 

K=16 

ISF=K 

REAL  COEFFICIENTS  --  SIGNAL  DEPENDENT  PARAMETERS: 

SZER  -  SIGMA  SUB  ZERO,  THE  SACLED  SOURCE  LEVEL 
(OR  SCALED  RAYLEIGH  DISTANCE) 

AZER  -  A  SUB  ZERO,  THE  ATTENUATION  OF  THE  FUNDAMENTAL  WITHIN 
THE  RAYLEIGH  DISTANCE 

B  -  WIDTH  OF  RADIAL  SAMPLE  SAPCE  RELATIVE  TO  RHO  SUB  ZERO 
NOTE:  THE  PROGRAM  USES  THE  SIGNAL  PARAMETERS  TO  DEDUCE 
RHO  SUB  ZERO. 

ALL  REFERENCES  TO  THE  'FUNDAMENTAL'  REFER  TO  THE  SIGNAL  WHICH 
EXECUTES  ONE  CYCLE  IN  THE  SAMPLE  WINDOW  AND  WHICH  HAS  AN  AMPLITUDE 
EQUAL  TO  THE  PEAK  AMPLITUDE  OF  THE  COMPOSITE  SIGNAL. 

—  MEDIUM  DEPENDENT  PARAMETERS: 

EMM  -  DISPERSION  PARAMETER 
TSUBC  -  RELAXATION  TIME 

~  PROCRAM  OPERATIONAL  PARAMETERS: 

SMAX  -  MAXIMUM  SCALED  RANGE  OF  INTEREST 
DSIG  -  DELTA  SIGMA,  THE  RANGE  INCREMENT 

DBLV  -  LEVEL  OF  THE  BEAM  EDGE  RELATIVE  TO  THE  CENTER  WHICH  WILL 
TRIGGER  RADIAL  EXAPNSION  OF  COORDINATES 
SIGCH  -  IF  GREATER  THAN  ZERO,  THE  RANGE  AT  WHICH  DELTA  SIGMA 
WILL  BE  CHANGED 

DSIGN  -  NEW  VALUE  FOR  DELTA  SIGMA  IF  A  CHANGE  IS  DESIRED 
EPRNT  -  RADIAL  DISTANCE(S)  OF  INTEREST  FOR  PRINTOUT 

READ( ICARD, 8020)  SZER,  AZER,  B 
READ! ICARD, 8020)  EMM,  TSUBC 

READ! ICARD, 8U20)  SMAX,  DSIG,  DBLV,  SIGCH,  DSIGN 
READ( ICARD, 8020)  (EPRNT(K) , K=1 ,8) 

CALL  L I NE NW (  ICASE,  IDAT  ) 

GAMMA  -  SCALED  ATTENUATION  COEFFICIENT. 

IRF  -  NUMBER  OF  RADIAN  FREQUENCY  COMPONENTS  IN  THE  DFT. 
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C  DE  -  NORMALIZED  RADIAL  DISTANCE  STEP,  DELTA  EPSILON. 

C  DEPS  -  DELTA  EPSILON,  USED  IN  THE  DIFFRA  SUBROUTINE.  THIS  PARAMETER 
C  FOLLOWS  THE  RADIAL  EXPANSION  OF  THE  COORDINATE  SYSTEM. 

C  STST  -  TEST  SIGMA  MAXIMUM;  DECREMENTED  TO  ALLOW  FOR  ROUNDOFF  ERROR 
C  SLEVEL  -  SPECTRUM  LEVEL  NORMALIZATION  FACTOR. 

C  TESTCH  -  SIGMA  VALUE  FOR  TEST  PURPOSES  AT  WHICH  DELTA  SIGMA  IS  TO 
C  BE  CHANGED.  IF  NO  CHANGE  IS  DESIRED,  IT  IS  TWICE  THE  MAXIMUM 

C  RANGE  OF  INTEREST,  WHICH  IS  NEVER  ATTAINED. 

C 

GAMMA  =  SZER/AZER 
IRF=INL/2 
JRF=IRF+1 
DE=B/FLOAT( IFB-1 ) 

DEPS=D£ 

STST=SMAX-0.5*DSIG 

SLEVEL=1. 

IF  (NCOMP  .NE.  0.)  SLEV EL=0 . 5 *SQ RT ( FL OAT ( I NL ) ) / FLO AT( NCOMP ) 
TESTCH=2 . *SMAX 

IF  (SIGCH  .GT.  0.)  TESTCH=SIGCH-0.5*DSIG 
C 

C  IFILE  -  DSRN  FOR  PERMANENT  RECORD 

C  NZP  -  NUMBER  OF  RANGE  STEPS  BETWEEN  OUTPUT  CYCLES. 

C  IFBO  -  NUMBER  OF  RADIAL  DISTANCE  STEPS  BETWEEN  PRINTOUT 
C  CYCLES 

C 

READ(ICARD,8000)  IFILE,  NZP,  IFBO 
MZP=0 
C 

C  TITLE  PAGE  INFORMATION. 

C 

WRITE(6,9010)  INL,  IFB,  ISF 
WRITE! 6,9170)  (IDIFFR(KK) ,KK=1,ISF) 

WRITEi 6,9020)  SZER,  AZER,  SMAX,  DSIG 
WRITE! 6,9030)  B,  DE,  GAMMA 
WRITE! 6,91 10)  EMM.  TSUBC 
WRITE(6,9330)  DBLV,  QLL.  RSPH 
IF  (NOLOSS)  WRITE! o,90A(j) 

IF  (NODISP)  WRITEl 6,9050) 

IF  1N0DIFF)  WRITE(6.9210) 

IF  (NOS PM)  WRITE(6,9220) 

IF  INONONL)  17RITEC 6,9260) 

IF  (SIGCH  .GT.  0.)  WRITE(6,9600)  SIGCH,  DSIGN 
600  CALL  LINENW(ICASE, IDAT) 

C 

C  SET  UP  THE  INITIAL  SIGNAL  WAVEFORMS. 

C 

CALL  UNIT2DC  U,  258,  40,  INL,  IFB,  B,  DMAG,  DE,  ICARD  ) 

C 

C  FILL  IN  THE  EPSILON  VALUE  ARRAY 
C 

FF=-DE 

DO  700  1=1, IFB 
FF=FF+DE 
EPSV( I)=FF 
700  CONTINUE 
C 

C  COMPUTE  THE  ATTENUATION  FACTORS 

C 

CALL  RALPHA(ATTF, INL, JRF, DSIG, GAMMA) 

C  NOTE:  PURGE  ATTF  ARRAY  OF  VALUES  LESS  THAN 

C  -  l.E-8;  THEY  WOULD  NOT  CONTRIBUTE  ANYWAY. 

C 

DO  800  1=2, JRF 

IF  (ATTF (I )  .LT.  l.E-8)  ATTF(I)=U. 

800  CONTINUE 
C 

C  COMPUTE  THE  DISPERSION  FACTORS. 

C 

CALL  CDISP(  SPMD,  JRF,  DSIG,  EMM,  TSUBC  ) 


C 
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C  CREATE  SPECTRAL  MODIFICATION  ARRAY  SPMD  IF  NECESSARY. 

C 

DO  1000  1=1, JRF 
CFC=1. 

IF  (NODISP)  GO  TO  900 

CFC=SPMD( I ) 

900  IF  (NOLOSS)  GO  TO  1000 
CFC=CFC*ATTF( I) 

1000  SPMD(I)=CFC 
C 

C  OUTPUT  THE  HEADER  INFORMATION  IF  REQUESTED  —  KREC  COUNTS  THE 
C  RECORDS. 

C 

CALL  SSWTCH(9,NORECD) 

NORECD= .NOT.NORECD 
IF  (NORECD)  GO  TO  1100 
KREC=1 

DF=DSIG*SZER 

WRITE(IFILE,8510)  KREC , INL, IFB, 1SF ,ICASE, IRF, NZP , IFBO , IDIFFR 
WRITE! IFILE, 8520)  I  DAT,  NODISP,  NOLOSS,  NODIFF 
WRITE(IFILE,8500)  SZER, AZER, SMAX, DSIG, B ,EMM,TSUBC ,DBLV , 

1  GAMMA,  DF,  ATTF,  SPMD 

C 

C  DISPLAY  INITIAL  FUNCTION  IF  REQUESTED. 

C 

1100  CALL  SSWTCH( 1 1 ,S 11 ) 

IF  ( .NOT . Si  1 )  GO  TO  1400 
CALL  SSWTCH(4 ,SW) 

IF  (.NOT.SW)  GO  TO  1300 
IDUM=3+(INL/10) 

DO  1200  1=1, IFB 

CALL  LINECK! IDUM, ICASE, IDAT) 

WRITE! 6, 9070)  I 

WRITE(6,9080)  (U(KK, I) , KK=1 , INL) 

1200  CONTINUE 

CALL  LINENW! ICASE, IDAT) 

1300  I DUM=3+ ( IRF / 5 ) 

CALL  LINECK! IDUM, ICASE, IDAT) 

WRITE! 6, 9120) 

WRITE! 6,9 130)  (SPMD(KK) , KK=1 , IRF) 

CALL  LINENW(ICASE, IDAT) 

C 

C  INITIALIZE  FOR  RANGE  STEP  PROCESSING. 

C  DS  -  PARAMETER  FOR  NONLINEAR  OPERATOR. 

C  NSTEP  -  ITERATION  COUNTER. 

C  SIGMA  -  SCALED  RANGE. 

C  BFAC  -  BEAM  WIDTH  NORMALIZATION  FACTOR. 

C  EFAC  -  ENERGY  NORMALIZATION  FACTOR. 

C 

1400  DS=FLOAT( INL )*DSIG/6. 2831853 
NSTEP-0 
SIGMA=0 . 

BFAC= 1 . 

EFAC=1 . 

FIF=FLOAT( INL) 

C 

C  INCREMENT  SIGMA,  THE  SCALED  RANGE. 

C 

1500  SIGMA=SICMA+DSIG 
C 

C  SEE  IF  DELTA  SIGMA  MUST  BE  CHANGED  --  IF  SO,  DO  SO. 

C 

IF  (SIGMA  .LT.  TESTCH)  GO  TO  1550 
RSIGCH=DSIGN/DSIG 
DO  1520  J= 1 , JRF 
SP I=SPMD( I ) 

SPMA=CABS( SP I ) **RSIGCH 

SPMP=RS ICCH*ATAN2 ( AI MAG( SP I ) , R  EAL( SP I ) ) 

SPMD( I ) =SPMA*CMPLX( COS( SPMP) ,SIN(SPMP) ) 

1520  CONTINUE 
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C 

C 

C 

1550 

C 

C 

C 

C 

C 

C 

C 


1600 

C 

C 

C 


1700 

C 

C 

C 

1300 

1900 

C 

C 

C 

2000 

C 

C 

C 


2100 

C 

C 

c 

2200 


2300 

2400 

C 

C 


TESTCH=2 .  *SMAX 

DSIG=DSIGN 

STST=SMAX-0.5*DSIG 

AFTER  (OR  WITHOUT)  CHANGE  OF  DELTA  SIGMA,  INCREMENT  COUNTERS. 

NSTEP=NSTEP+1 

MZP=M2P+1 

APPLY  THE  NONLINEAR  OPERATOR  WITHIN  THIS  LOOP  TO  EACH  SET  OF 
WAVEFORMS  SAMPLES  AT  THE  RADIAL  DISTANCE  INCREMENT  COUNTED  BY  J. 

DO  2400  J=1 , IFB 

MOVE  THE  J-TH  RADIAL  WAVEFORM  SAMPLE  INTO  UR  ARRAY. 

DO  1600  K=1 , INL 
UR(K)=U(K,J) 

CONTINUE 

IF  (NONONL)  GO  TO  1800 
NON-LINEAR  OPERATOR  ON  UR  INTO  UZ. 


DX=0. 

DO  1700  11=1, INL 
DX=DX+1 . 

DL=DS*UR( II) 

XT=DX+DL 

IF  (XT  .LT.  0.)  XT=XT+FIF 

M0=XT 

Ml=MO+l 

DL=XT-FLOAT(MO) 

IF  (MO  .LT.  1)  M0=M0+INL 
IF  (MO  .GT.  INL)  M0=M0-INL 
IF  (Ml  .LT.  1)  M1=M1+INL 
IF  (Ml  .GT.  INL)  M1=M1- INL 
UZ(II)=UR(M0)*(1.-DL)  +  UR(Ml)*DL 
CONTINUE 
GO  TO  2000 

BYPASS  NONLINEAR  OPERATOR. 

DO  1900  11=1 .INL 
UZ(II)=CMPLX(UR(II),0.) 

CONTINUE 

DIRECT  FOURIER  TRANSFORM  —  TIME  TO  FREQUENCY  DOMAIN. 

CALL  DFT  (  UZ,  W,  ISQ,  INL,  2  ) 

APPLY  ATTENUATION  (DISPERSION). 

IF  (NOS PM)  GO  TO  2200 
DO  2100  K=2. IRF 
UZ(K)=UZ(K)*SPMD(K) 

CONTINUE 

UZ( JRF )=UZ( JRF )*CABS(SPMD( JRF ) ) 

RESTORE  RADIAN  FREQUENCY  SAMPLES  TO  U  ARRAY. 

MM=0 

DO  2300  K=l, JRF 
MM=MM+1 

U(MM, J)=REAL(UZ(K) ) 

MM=MM+1 

U(MM, J )=AIMAG(UZ(K) ) 

CONTINUE 

CONTINUE 

SEE  IF  DIFFRACTION  IS  TO  BE  BYPASSED. 
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C 

IF  (NODIFF)  GO  TO  3100 
C 

C  PERFORM  THE  ZERO-FREQUENCY  'DIFFRACTION' 

C  FIRST  PICK  UP  THE  ZERO-FREQUENCY  (DFT  CELL  NUMBER  1)  COMPONENTS; 
C  THEN  'DIFFRACT'  THEM  USING  SUBROUTINE  DIFFRO. 

C 

DO  2500  1=1, IFB 
UO(I)=CMPLX(U(l,I),U(2,I)) 

2500  CONTINUE 

CALL  DIFFR0(U0 , IFB , SIGMA, DS IG ,SZER) 

DO  2550  1=1, IFB 
U(l,I)=REAL(UO(I)) 

U(2.I)=AIMAG(UO(I>) 

2550  CONTINUE 
C 

C  PERFORM  THE  DIFFRACTION  CALCULATION 

C  FOR  EACH  OF  THE  RADIAN  FREQUENCIES,  SAMPLES  AT  EACH  RADIAL 

C  DISTANCE  ARE  MOVED  INTO  ARRAY  UO. 

C 

FI=0 . 

LL=1 

MM=2 

DO  3000  1=2, JRF 
LL=LL+2 
MM=M>l+2 
FI=FI+1 . 

C 

C  SELECT  THE  RADIAL  SAMPLES  AT  FREQUENCY  I,  MOVE  INTO  UO  ARRAY. 

C 

DO  2600  J= 1 , IFB 
U0(J)=CMPLX(U(LL,J),U(MM,J)) 

2600  CONTINUE 
C 

C  PERFORM  THE  DIFFRACTION. 

C 

CALL  DIFFRA  (  UO ,  IFB,  DSIG,  DEPS,  SZER,  FI,  QLL,  SIGMA,  RSP11  ) 
C 

C  RESTORE  THE  DIFFRACTED  SPECTRAL  SAMPLES  TO  THE  U  ARRAY. 

C 

DO  2700  J=1 , IFB 
U(LL , J )=REAL( UO ( J  )  ) 

U(MM, J)=AIMAG(UO(J)) 

2700  CONTINUE 
3000  CONTINUE 
C 

C  OUTPUT  THE  SPECTRUM  IF  REQUESTED. 

C 

3100  IF  (MZP  ,NE.  NZP )  GO  TO  3900 
IF  (NORECD)  GO  TO  3300 
KREC=KREC+1 

WRITE(IFILE,8530)  KREC, ICASE, DEPS , SIGMA 
DO  3200  JJ=  J. ,  IFB 

WRITE  (IFILE.8540)  (U(KK, JJ) , KK= 1 , INL) 

3200  CONTINUE 
C 

C  DETERMINE  THE  AXIAL  LEVELS  AND  THE  BEAM  WIDTHS. 

C 

3300  LL=-1 

MM=0 

SST=SIGMA*0. 2399/(SZER*DEPS) 

Frs= 1 m 

DO  3340  1=1, JRF 

LL=LL+2 

MM=MM+2 

IF  (I  .GT.  1)  FRS=FLOAT(I— 1 ) 

NAXIS=1+IF IX(SST/FRS) 

DO  3320  J=1 , IFB 

DMAG( J )=SQRT(U(LL, J)**2+U(MM,J )**2 )/ SLEVEL 
3320  CONTINUE 
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FMAX=DMAG( 1 ) 

IF  (NAX1S  .LT.  2)  GO  TO  3330 
CALL  MAXFCN(DMAG, NAXIS, FMAX, INDEX) 

3330  AXISLV(I )=FMAX 
3340  CONTINUE 
LL=-1 
MM=0 

DO  3600  1=1, JRF 

LL=LL+2 

MM=MM+2 

DO  3400  J= 1 , IFB 

DMAG( J )=SQRT(U(LL, J )**2+U(MM, J )**2 )/ SLEVEL 
3400  CONTINUE 

DO  3500  J=1 .3 
TEST=AXISLV(I)*ROOTV(j) 

CALL  FINDl ( 1 , IFB, EPSV.DMAG, TEST, POINT, IDUM) 

BW( J ,I)=0. 

IF  (IDUM  .EQ.  I)  BW(J,I)=POINT*BFAC 
3500  CONTINUE 
3600  CONTINUE 

CALL  LINECK( 2, ICASE, IDAT) 

WRITE(6,9300)  BFAC 
C 

C  BYPASS/FINISH  SPECTRUM  OUTPUT. 

C 

3900  KFB0=1 

JDUM=2+( (IRF+9)/ 10) 

C 

C  SEE  IF  RADIAL  COORDINATE  EXPANSION  IS  REQUIRED.  EACH  OF  THE 

C  NON-ZERO  FREQUENCIES  IN  ARRAY  IDIFFR  IS  TESTED  TO  SEE  IF  THE  BEAM 

C  HAS  BROADENED  SUFFICIENTLY  TO  REQUIRE  RADIAL  EXPANSION. 

C 

IF  (( SIGMA/ SZER)  .LT.  RDIFF )  GO  TO  4500 
SQM=-1. 

DO  4100  1=1, ISF 
IDI=2*IDIFFR(I)+1 

SQS=SQRT( (U( IDI . IFB) **2+U( IDI+1 , IFB ) **2 )/ 

1  (U( IDI , 1 )**2+u{ IDI+1 , 1 )**2 ) ) 

IF  (SQS  .GT.  SQM)  SQM=SQS 
4100  CONTINUE 

IF  (SQM  -LT.  DBLV)  GO  TO  4500 
C 

C  EXPAND  RADIAL  COORDINATES. 

C 

LL=~1 

MM=0 

DO  4400  J=1 , JRF 

LL=LL+2 

MM=MM+2 

DO  4200  K=1 , IFB 
UO(K)=CMPLX(U(LL,K) ,U(MM,K) ) 

4200  CONTINUE 

CALL  EXPAND(U0 , IFB.XFAC) 

DO  4300  K= 1 , IFB 
U(LL, K)=REAL(U0 (K ) ) 
u(mm,k)=aimag(uo(k)) 

4300  CONTINUE 
4400  CONTINUE 
C 

C  ADJUST  BFAC,  THE  BEAM  WIDTH  NORMALIZATION,  DEPS,  THE  RADIAL 
C  DISTANCE  INCREMENT,  AND  EFAC,  THE  BEAM  ENERGY  NORMALIZATION 
C  FACTOR,  TO  REFLECT  THE  EXPANDED  RADIAL  COORDINATES. 

BFAC=BFAC*XFAC 
DEPS=DEPS*XFAC 
EFAC=EFAC*SQRT ( XFAC) 

IDUM=2 

CALL  LINECK( IDUM, ICASE, IDAT) 

WRITE(6,9290)  SIGMA,  BFAC,  EFAC 


C 


C  SELECT  EACH  OF  THE  IFB  RADIAL  SAMPLE  SPECTRA  IN  SUCCESSION 
C  AT  SPATIAL  FREQUENCY  J. 

C 

4500  DO  5100  J=1 , IFB 
C 

C  SELECT  RADIAN  FREQUENCY  SAMPLES  AT  FREQUENCY  K. 

C 

LL=-1 

MM=0 

DO  4600  K=1,JRF 

LL=LL+2 

MM=MM+2 

U0(K)=CMPLX(U(LL, J ) ,U(MM, J ) ) 

4600  CONTINUE 
NE=JRF 
NB=JRF+1 

DO  4650  K=NB, INL 
NE=NE— l 

UO(K)=CONJG(UO (NE) ) 

4650  CONTINUE 

IF  (J  .NE.  KFBO)  GO  TO  4800 
IF  (MZP  .NE.  NZP)  GO  TO  4800 
C 

C  PICK  UP  THE  MAGNITUDE 
C 

DO  4700  KK=1 , IRF 
DMAG(KK)=CABS(UO (KK) )/SLEVEL 
4700  CONTINUE 
C 

C  DISPLAY  SPECTRUM  IF  REQUESTED. 

C 

KFBO=KFBO+IFBO 

CALL  LINECK( JDUM, ICASE, IDAT) 

WRITE(6,9190)  J,  NSTEP 
WRITE(6,9080)  (DMAG(KK) , KK=1 , IRF) 

C 

C  INVERSE  FOURIER  TRANSFORM  THE  DATA  AT  EACH  RADIAL  SAMPLE  DISTANCE. 
C  THIS  SECTION  CONVERTS  THE  NON-NEGATIVE  FREQUENCY  COMPONENTS 
C  STORED  IN  'U'  INTO  REAL  WAVEFORMS  AT  EACH  OF  THE  RADIAL  DISTANCES. 
C 

4800  CALL  DFT(U0 ,  W,  ISQ,  INL,  -2  ) 

C 

C  RESTORE  THE  WAVEFORMS  INTO  THE  U  ARRAY. 

C 

4900  DO  5000  K=1 , INL 

U(K, J )=REAL(U0 (K) ) 

5000  CONTINUE 
5100  CONTINUE 
C 

C  OUTPUT  THE  AXIAL  LEVELS  AND  THE  BEAM  WIDTHS. 

C 

IF  (MZP  .NE.  NZP)  GO  TO  5200 

I D  UM=4+ ( J  RF+9 ) / 1 0 

CALL  LINECK(IDUM, ICASE. IDAT) 

WR IT E( 6,9610)  NSTEP,  SIGMA 
WRITE(6,9080)  (AXISLV(KK) , KK=1 , JRF ) 

IDUM=2+(  J  RF+9)/  J.0 

CALL  LINECK(IDUM, ICASE, IDAT) 

WRITE(6,9140) 

WRITE(  6 , 9080)  ( BW(  1. , KK)  ,  KK=  1 ,  JRF ) 

CALL  LINECKdDUM,  ICASE,  IDAT) 

WKITE( 6,9150) 

WRITER 6 , 9080)  (BW( 2 ,KK) , KK= 1 , JRF ) 

CALL  LINECKdDUM,  ICASE,  IDAT) 

WRITE(6,9160) 

WRITE(6,9080)  (BW(3,KK) ,KK=1 ,JRF) 

C 

C  RETURN  FOR  NEXT  RANGE  STEP. 

C 

5200  R=SICMA/SZER 
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C 

C 

C 

7000 


7040 


7080 

7120 


7160 


7200 

C 

C 

C 

8000 

8020 

C 

C 

C 

8500 

8510 

8520 

8530 

8540 

C 

C 

C 

9010 

9020 

9030 

9040 

9050 

9060 

9070 

9080 

9090 

9110 

9120 

9130 

9140 

9150 


IF  (SIGMA  .GE.  ST ST)  GO  TO  7000 
IF  (MZP  .GE.  NZP)  MZP=0 
CALL  LINECK( 1 , I CASE, IDAT) 

WRITE(6,906Q)  NSTEP,  SIGMA,  R 
GO  TO  1500 

PRINT  WAVEFORM(S)  IF  REQUESTED;  SEE  IF  TASK  IS  FINISHED. 

CALL  LINENW  (  ICASE,  IDAT  ) 

BB=BFAC*B 
IDUM=(INL/10)+4 
DO  7080  1=1,8 
EPI=EPRNT(I) 

IF  (EPI  .LT.  0.)  CO  TO  7120 
IF  (EPI  .GT.  BB)  GO  TO  7120 
E0=1.+(£PI/D£PS) 

LE=IFIX(EO) 

FEO=EO-FLOAT(LE) 

GE0=1.-FE0 
ME=LE+1 

DO  7040  J=1 , INL 

DMAG( J )=U( J ,LE) *GEO+U( J , ME ) *FEO 
CONTINUE 
CALL  LINECK 
WRITE(6,9240 
WRITE(6,9080 
WRITE( 6, 9080 
CONTINUE 
READ(ICARD,8020)  SMAX 
IF  (SMAX  .LT.  0.)  GO  TO  7160 
READ! ICARD, 8020)  (EPRNT(K) , K=1 ,8) 

STST=SMAX-0.5*DSIG 
GO  TO  1500 

IF  (NORECD)  GO  TO  100 
KREC=-1 

WRITE(IFILE,8530)  KREC, ICASE, DEPS, SIGMA 

GO  TO  100 

WRITE(6,9090) 

STOP 

INPUT  FORMATS. 

FORMAT! 1615) 

FORMAT( 8F 10. 1 ) 

OUTPUT  FORMATS  FOR  DATA  STORAGE  FILE. 

FORMAT! 1P8E10.3) 

FORMAT! 1018) 

FO  RMAT{ 2A4 , 3L 1 ) 

FORMAT! 2110 , 1P6E10 . 3) 

FORMAT! 1P32E10 .3) 

OUTPUT  FORMATS  FOR  PRINTED  RECORD. 

FORMAT! 'OINL' ,14, '  IFB',13,'  ISF',13) 

FORMAT!'0  SZER'. 1PE12.3,'  AZER' ,E12.3, '  SMAX',E12.3, 

1  '  DSIC',Eli.3) 

FORMAT! '0  B', 1PEI2.3,'  DE',E12.3,'  GAMMA' ,E12. 3) 

FORMAT! '0 LINEAR  LOSSES  SUPPRESSED.') 

FORMAT! 'ODISPERSION  SUPPRESSED.') 

FORMAT!'  STEP', 15,'  SIGMA' , 1PE12.3.'  R',EJ2.3) 

FORMAT!'0 INITIAL  FUNCTION  AT  RADIAL  STEP', 14) 

FORMAT!'  ' , 1P1 0E12 . 3) 

FORMAT! 'OPROGRAM  REACHED  NORMAL  TERMINATION.') 

FORMATrO  EMM'  1PE12.3, '  TSUBC' ,E12. 3) 

FORMATTO  SPECTRUM  MODIFICATION  ARRAY  (ATTENUATION/DISPERSION):') 
FORMAT!'  '  1P10E13.4) 

FORMAT!'0  3  DB  BEAM  WIDTHS:') 

FORMAT! '0  6  DB  BEAM  WIDTHS:') 


(  IDUM,  ICASE,  IDAT  ) 

]  NSTEP,  SIGMA,  R,  EPI 

)  (DMAG(K) , K=1 ,INL) 
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9160  FORMAT! '010  DB  BEAM  WIDTHS:') 

9170  FORMAT! 'OTEST  FOR  DIFFRACTION: 1615) 

9190  FORMAT! '  0  SPECTRUM  SAMPLE'  14/  AFTER  OPERATOR  STEP', 15) 

9210  FORMAT! 'ODIPFRACT ION  OPERATOR  SUPPRESSED.') 

9220  FORMAT! ' ORADIAN  FREQUENCY  SPECTRUM  MODIFICATION  SUPPRESSED.') 

9240  FORMAT! 'OSTEP' ,15,'  SIGMA,  R: ' , 1P2E12 .3 , ' ;  RADIUS:', 

1  E12.3) 

9260  FORMAT! 'ONONLINEAR  OPERATOR  SUPPRESSED.') 

9290  FORMAT! 'OAT  SIGMA=' , 1PE9.2,'  RADII  EXPANDED  TO',E9.2,'  AND  BEAM  EN 
1ERGY  TO'  E9.2) 

9300  FORMAT! 'OCURRENT  BFAC' , 1PE12 .4) 

9330  FORMAT! '0  DBLV'  1PE12 .3/  QLL',E12.3,'  RSPH',E12.3) 

9600  FORMAT! 'OAT  SIGMA=' , 1PE12 .4, '  DELTA  SIGMA  BECOMES' ,E12 .4) 

9610  FORMAT! 'OAXIAL  BEAM  LEVELS  AT  STEP' ,15,'  SIGMA' , 1PE12.4) 

END 


SUBROUTINE  DIFFRO  (  UO ,  N,  SIGMA,  DSIG,  SZER  ) 

C 

C  THIS  SUBROUTINE  PERFORMS  THE  'ZERO-FREQUENCY'  DIFFRACTION  FOR 
C  THE  PRPGTS  NONLINEAR  WAVE  EQUATION  SOLUTION  PROGRAM. 

C 

C  RATIONALE: 

C  SPHERICAL  OR  1/R  SPREADING  IS  IMPOSED.  THE  DIFFUSION-IN-FREQUENCY 

C  ASSUMES  THE  FORM  OF  LAPLACE'S  EQUATION  !DEL  SQUARED  PHI  =  0)  IN  THE 

C  LIMIT  AS  FREQUENCY  GOES  TO  ZERO.  A  SOLUTION  TO  LAPLACE'S  EQUATION 
C  IS  1/R. 

C 

C  DEFINITIONS  OF  ARGUMENTS 

C  UO  -  COMPLEX  DATA  ARRAY  TO  BE  'DIFFRACTED' 

C  N  -  NUMBER  OF  POINTS  IN  UO 

C  SIGMA  -  SCALED  RANGE  ATTAINED 

C  DSIG  -  SCALED  STEP  SIZE 

C  SZER  -  SCALED  DIFFRACTION  PARAMETER  !RAYLEIGH  DISTANCE) 

C 

C  SOLI  DEO  GLORIA  5>  XI  1979  —  F  S  MCKENDREE 

C 

COMPLEX  F,  UO!N) 

COMPLEX  UQ(50) 

C 

C  RN  -  NEW  RANGE 

C  RO  -  OLD  RANGE 

C  F  -  RANGE  DECREMENT  FACTOR 

C 

RN=SIGMA/SZER 

R0=! SIGMA- DS IG)/ SZER 

F=CSQRT!CMPLX! 1 . ,-RO)/CMPLX< I . ,-RN) ) 

C 

C  APPLY  RANGE  DECREMENT  FACTOR  TO  THE  DATA 
C 

DO  1000  1=1 ,N 
UQ!l)=UO!l)*F 
1000  CONTINUE 
C 

C  SPREAD  THE  DATA  SPHERICALLY 
C 

DG=! 1 .+SIGMA )/ ( 1 ,+SIGMA+DS IG) 

EG=0. 

UO(i)=UQ(l) 

DO  2000  1=2, N 
EG=EG+DG 
J=IFIX(EG)+1 
K=J+1 

FC=EG-FL0AT!J-1) 

HG=J .-FG 

UO(I)=UQ(J)*HC+UQ(K)*FG 
2000  CONTINUE 
RETURN 
END 
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SUBROUTINE  EXPAND  (  P,  N,  XFAC  ) 

C 

C  100  PERCENT  RADIAL  EXPANSION  —  INVOLVES  NO  INTERPOLATION. 
C 

COMPLEX  P(N) ,  CF/ (0 .  ,0 . )/ 

M-l 

K=1 

1000  M=M+2 

IF  (M  .GT.  N)  GO  TO  2000 

K=K+1 

P(K)=P(M) 

GO  TO  1000 
2000  K=K+1 

IF  (K  .GT.  N)  GO  TO  8888 
DO  3000  I=K,N 
P(I)=CF 
3000  CONTINUE 
8888  XFAC=2 . 

RETURN 

END 

SUBROUTINE  MAXFCN  (  F,  N,  FMAX,  INDEX  ) 

C 

C  F  -  REAL  ARRAY  TO  BE  SEARCH  FOR  MAXIMUM 
C  N  -  LARGEST  SAMPLE  NUMBER  TO  BE  SEARCHED 
C  FMAX  -  RETURNED  AS  THE  LARGEST  VALUE  IN  F 
C  INDEX  -  RETURNED  S  THE  LOCATION  OF  FMAX  IN  F 
C 

C  SOLI  DEO  GLORIA  8  III  1981  —  F  S  MCKENDREE 

C 

REAL  F(N) 

FMAX  =  F(l) 

INDEX=1 

IF  (N  .LT.  2)  GO  TO  8888 
DO  1000  1=2, N 

IF  (F(I)  .LE.  FMAX)  GO  TO  1000 
FMAX=F(I) 

INDEX=I 

1000  CONTINUE 
8888  RETURN 
END 


SUBROUTINE  UNIT2D  (  U,  IN,  IF,  INL,  IFB,  RHOZER,  RHO,  DR,  ICARD  ) 

THIS  SUBROUTINE  PROVIDES  INITIALIZATION  OF  THE  DATA  ARRAY  IN  THE 
TIME  DOMAIN  AS  REQUIRED  BY  PRPGT8.  THREE  BEAM  TYPES  ARES  PROVIDED, 
INCLUDING  THOSE  USED  BY  BAKHVALOV  ET.  AL. 

SOLI  DEO  GLORIA  21  DECEMBER  1979  —  F  S  MCKENDREE 

DEFINITIONS  OF  ARGUMENTS: 

U  -  REAL  ARRAY  OF  TWO  DIMENSIONS  IN  BY  IF,  WHICH  WILL  CONTAIN 
THE  ORIGINAL  TIME  WAVEFORMS  UPON  RETURN. 

IN. IF  -  DIMENSIONS  OF  THE  U  ARRAY. 

INL  -  NUMBER  OF  SAMPLES  OF  THE  WAVEFORM  AT  A  GIVEN  RADIAL 
DISTANCE  —  I.  E.,  THE  NUMBER  OF  SAMPLE  ROWS  ALONG  THE 
DIRECTION  OF  PROPAGATION. 

IFB  -  NUMBER  OF  SAMPLES  ACROSS  THE  WIDTH  OF  THE  BEAM  AT  EACH 

TIME  SAMPLE  —  I.  E.,  THE  NUMBER  OF  SAMPLE  COLUMNS  ACROSS 
THE  RADIAL  DIRECTION. 

RHOZER  -  RHO  SUB  ZERO,  NO  LONGER  USED  BUT  INCLUDED  FOR  PURPOSES 
OF  COMPATABlLlTY  WITH  EARLIER  ROUTINES. 

RHO  -  REAL  ARRAY  USED  FOR  TEMPORARY  STORAGE  OF  THE  BEAM  PROFILE 
AMPLITUDE. 

DR  -  RADIAL  DISTANCE  STEP  SIZE. 

ICARD  -  DSRN  TO  BE  USED  FOR  CONTROL  INPUT  TO  'UNIT2D' . 

REAL  RHO(IF),  U(IN,IF) 
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C 

C  ZERO  THE  U  ARRAY. 

C 

DO  500  1=1, IF 
DO  500  J=1 , IN 
U(J,I)=0. 

500  CONTINUE 
C 

C  READ  A  DATA  'CARD'.  THIS  CARD  SPECIFIES  THE  FREQUENCY,  AMPLITUDE, 

C  AND  PHASE  OF  UP  TO  TWO  COMPONENTS. 

C  FI  LESS  THAN  ZERO  CAUSES  THE  SUBROUTINE  TO  NORMALIZE  THE 

C  SAMPLE  SET  TO  UNIT  PEAK  AMPLITUDE  AND  RETURN. 

C  IF  THE  PHASE  IS  IN  THE  RANGE  0  TO  360  DEGREES,  A  GAUSSIAN  BEAM 

C  PROFILE  IS  USED;  IF  IN  THE  RANGE  360  TO  720  DEGREES,  A  BOXCAR  PROFI 

C  LE  IS  USED;  AND  IF  IN  THE  RANGE  720  TO  1080  DEGREES,  A  FOURTH  ORDER 

C  POLYNOMIAL  IS  USED. 

C 

WRITE(6, 9100)  ICARD 
WRITE (6, 9000) 

1000  READ (ICARD, 80 00)  FI,  Al .  PI,  F2 ,  A2 ,  P2 
IF  (FI  .LT.  0.)  GO  TO  4000 
R=-DR 

IF  (PI  .LT.  360.)  GO  TO  1200 
IF  (PI  .LT.  720.)  GO  TO  1800 
IF  (PI  .LT.  1080.)  GO  TO  2200 
C 

C  GAUSSIAN  BEAM  —  DEFAULT  OR  PHASE  <  360 
C 

1200  WRITE(6,9020) 

DO  1500  1=1, IFB 
R=R+DR 

RHO(l)=EXP(-(R*R) ) 

1500  CONTINUE 

GO  TO  3000 
C 

C  UNIFORM  EXCITATION  TO  R=1  —  PHASE  <  720 
C 

1800  Pl=Pl-360. 

WRITE(6,9040) 

DO  2000  1=1, IFB 
R=R+DR 
RHO(I)=l . 

IF  (R  .GT.  1.)  RHO(I )=0 . 

2000  CONTINUE 

GO  TO  3000 
C 

C  FOURTH  ORDER  POLYNOMIAL  BEAM 
C 

2200  Pl=Pl-720. 

WRITE(6,9060) 

DO  2500  1=1, IFB 
R=R+DR 

RHO(I)=l./((l.+R*R)**2) 

2500  CONTINUE 
C 

C  HAVING  THE  BEAM  PROFILE  IN  R,  INSERT  THE  FREQUENCY  COMPONENTS. 

C 

3000  TAU=6 .28  31 8 53/ FLOAT ( INL) 

WRITE(6,9080)  FI,  Al,  Pi,  F2,  A2,  P2 
P1=P1*0. 01745329 
P2=P2*0. 01745329 
T=-TAU 

DO  3500  1=1, INL 
T=T+TAU 

S=Al *COS ( F 1 *T+P1 )+A2*COS ( F2  *T+P2 ) 

DO  3200  J= 1 , IFB 
U(I,J)=U(I, J)+S*RHO(J) 

3200  CONTINUE 
3500  CONTINUE 

IF  (F2  .LT.  0.)  GO  TO  4000 
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GO  TO  1000 
C 

C  NORMALIZE  TO  UNIT  PEAK  AMPLITUDE  —  FIRST  FIND  THE  MAXIMUM,  THEN 
C  ADJUST  THE  AMPLITUDES  ACCORDINGLY. 

C 

4000  AMAX=-1. 

DO  4400  1=1, INL 
DO  4200  J=1 , IFB 
A=ABS(U(I,J)) 

IF  (A  .GT.  AMAX)  AMAX=A 
4200  CONTINUE 
4400  CONTINUE 

AMAX=I./AMAX 
DO  4800  1=1, INL 
DO  4600  J=l, IFB 
U(I.J)=U(I,J)*AMAX 
4600  CONTINUE 
4800  CONTINUE 
RETURN 

8000  FORMAT(8F 10.1) 

9000  FORMAT! ' 0UNIT2D  Fl' , 10X, 'A1 ' , I0X , ' PI ' , 10X, 'F2 ' , 10X, 'A2 ' , 10X , 

9020  FORMAT! 'OGAUSS IAN') 

9040  FORMAT! 'OUNIFORM' ) 

9060  FORMAT! '04TH  ORDER') 

9080  FORMAT! '0' , 10X, 1P6E12 . 3) 

9100  FORMAT! 'ODATA  TAKEN  FROM  FILE ',15) 

END 


SUBROUTINE  FINDl  !  NB,  NE ,  ORDI ,  ABSC,  ROOT,  POINT,  ITEST  ) 

C  'ROOT'  FINDING  SUBROUTINE.  THE  ARRAY  ABSC  IS  SEARCHED  BETWEEN 
C  ELEMENTS  NB  AND  NE  TO  FIND  THE  VALUE  ROOT.  POINT  IS  RETURNED 
C  AS  THE  LINEARLY-INTERPOLATED  VALUE,  DETERMINED  FROM  ARRAY  ORDI, 

C  AT  WHICH  THE  FUNCTION  IN  ABSC  WILL  ATTAIN  THE  VALUE  ROOT.  ITEST 

C  IS  RETURNED  AS  ZERO  IF  NO  SUCH  VALUE  IS  FOUND  AND  AS  I  IF  SUCH 

C  A  VALUE  IS  FOUND. 

C 

C  DESIGNED/CODED  BY  F  S  MCKENDREE,  31  V  1978. 

C 

REAL  ORDl!NE),  ABSC!NE) 

C 

C  DETERMINE  IF  'ROOT'  IS  ATTAINED. 

C 

K=NB 

aold=absc!k) 

1000  K=K+1 

anew=absc!k) 

IF  (AOLD  .le.  root  .and.  anew  .ge.  ROOT)  GO  TO  2000 
IF  (AOLD  .GE.  ROOT  .AND.  ANEW  .LE.  ROOT)  GO  TO  2000 
AOLD=ANEW 

IF  (K  .LT.  NE)  GO  TO  1000 
POINT=0. 

ITEST=0 
8888  RETURN 
C 

C  FIX  VALUE  OF  'POINT'  CORRESPONDING  TO  'ROOT'. 

C 

2000  ITEST=1 

IF  (ANEW  ,EQ.  AOLD)  GO  TO  3000 
ORDO=ORDI( K-l ) 

ORDN=ORDl(K) 

POINT  =  ORDO  +  (((ROOT-AOLD)/(ANEW-AOLD))*(ORDN-ORDO)) 

GO  TO  8888 
C 

C  EQUAL  AMPLITUDES  BOTH  EQUAL  TO  ROOT. 

C  CHOOSE  A  VALUE  FOR  POINT  MIDWAY  IN  THE  CORRESPONDING  ORDINATES. 

C 

3000  POINT=(ORDI(K- I )+ORDI (K) )*0. 5 
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GO  TO  8888 
END 


SUBROUTINE  LINECK  (  LINES,  ICASE,  IDAT  ) 

CLEAN  PRINTER  CONTROLLER. 

REVISED  9  OCTOBER  1978  TO  INCLUDE  TITLE  ARRAY. 

LINECK(LINES, ICASE, IDAT)  -  DETERMINES  IF  'LINES'  MAY  BE  ADDED 
TO  THE  CURRENT  PAGE;  IF  NOT,  INCREMENT  PAGE  COUNTER  AND  MOVE  TO 
HEAD  OF  NEXT  PAGE,  WRITING  THE  TITLE  ARRAY,  'ICASE'  AND  THE  DATE 
FIELD  'IDAT'. 

LINETT( DUMMY)  -  MOVE  THE  TITLE  ARRAY  POINTED  TO  BY  'DUMMY'  INTO 
THE  'TITLE'  ARRAY. 

LINENW(ICASE, IDAT)  -  PAGE  IMMEDIATELY,  DISPLAYING  'ICASE'  AND 
'IDAT';  LINE  COUNTER  'LOCAL'  BECOMES  2. 

SOLI  DEO  GLORIA  5  OCTOBER  1978. 


INTEGER  NPAGE/0/ ,  LOCAL/60/,  IDAT(2) 
INTEGER  TITLE( 16 ) ,  DUMMY(16) 


SEE  IF  PAGING  IS  REQUIRED  --  IF  SO, 

LOCAL=LOCAL+LINES 
IF  (LOCAL  .LT.  60)  GO  TO  1000 
LOCAL=LINES+2 
700  NPAGE=NPAGE+1 

WRITE(6,9000)  TITLE,  ICASE,  IDAT,  NPAGE 
1000  RETURN 


ENTRY  LINETT (DUMMY) 

DO  2000  1=1,16 
TITLE( I )=DUMMY ( I ) 

2000  CONTINUE 
GO  TO  1000 

ENTRY  LINENW  (  ICASE, 

LOCAL=2 

GO  TO  700 


SET  UP  THE  TITLE  FIELD 

PAGE  IMMEDIATELY 
IDAT  ) 

THE  TITLE'S  FORMAT 


PAGE 


9000  FORMAT('l',16A4,'  CASE', 18,'  DATE  ',2A4,'  PAGE', I  4) 
END 


SUBROUTINE  RALPHA(ATTF,INL, IRF.DSIGMA, GAMMA) 

COMPUTE  THE  ATTENUATION  FACTORS  FOR  SPECIFIED  INPUT. 

REAL  ATTF(IRF) 

ATTF( 1 )=1 . 

K=0 

DSG= DSIGMA/GAMMA 
DO  1000  1=2, IRF 
K=K+1 

E=-FLOAT(K*K)*DSG 
IF  (E  .LT.  -60.)  E=-60 . 

ATTF(I)=EXP(E) 

1000  CONTINUE 
RETURN 
END 


SUBROUTINE  CDISP  (  DSP,  IRF,  DSIG,  EMM,  TSUBC  ) 
EVALUATE  THE  DISPERSION  COEFFICIENT. 
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DSP  -  ARRAY  TO  STORE  THE  COEFFICIENTS. 
INL  -  NUMBER  OF  FOURIER  COMPONENTS. 

DSIG  -  RANGE  STEP  PARAMETER. 

EMM  -  DISPERSION  PARAMETER. 

TSUBC  -  RELAXATION  PARAMETER. 

SOLI  DEO  GLORIA  OCTOBER  MCMLXXVIII. 

COMPLEX  DSP(IRF),  CA 

DSP( 1 )=1 . 

FI=0. 

FT=0 . 

DM=-DS IG*EMM 
DO  1000  1=2, IRF 
FI=FI+1 . 

FT=FT+TSUBC 

CA=DM* ( F  I**2)/CMPLX(1..FT) 

IF  (CABS(CA)  ,LT.  l.E-6)  GO  TO  800 
DSP(I)  =  CEXP(CA) 

GO  TO  1000 
DSP(I)=1. 

1000  CONTINUE 

DSP ( IRF )=REAL(DSP( IRF ) ) 

RETURN 
END 


SUBROUTINE  DIFFRA(W,N,DS,DEPS,SZER,FI,QLL,SIGMA,RSPH) 

THIS  SUBROUTINE  PERFORMS  THE  DIFFUSION-IN-FREQUENCY  CALCULATIONS 
BY  MEANS  OF  A  GAUSS  ELIMINATION  SOLUTION  OF  THE  IMPLICIT  CRANK- 
NICHOLSON  FORMULATION  OF  THE  PROBLEM. 

SOLI  DEO  GLORIA  29  OCTOBER  1980  —  F  S  MCKENDREE 

ADDITION  6  NOVEMBER  1980  —  CODE  FOR  SPHERICAL  SPREADING  ADDED 
AFTER  THE  DESIGN  OF  THIS  DATE. 

COMPLEX  W(N),  P(50) ,  AJ ( 50 ) ,  BJ(50),  CJ(50),  AF( ;0) ,  SF(50),  D(50) 

COMPLEX  Al,  A2,  Bl,  B2,  Rl,  A3,  A4,  A5 

COMPLEX  CG 

REAL  B/l ./  ,  C/2./ 

REAL  SSI/ 1 .065/ ,  SS2/7.000/,  SS3/1.750/,  SS4/2.000/ 

TEST  RANGE  FOR  SPHERICAL  SPREADING. 

RG=SIGMA/(SZER*FI) 

IF  (RG  .GT.  RSPH)  GO  TO  5000 

SET  UP  NUMBER  OF  ITERATIONS. 

T=C*DS*QLL/(0.5*FI*SZER*(DEPS**2)> 

NT=IFIX(T) 

IF  (NT  .LT.  1)  NT=1 
DSIG=0.25*DS/FL0AT(NT) 

SET  UP  FOR  ITERATION. 

Al=CMPLX ( 0. , DS IG/ ( 2 . *FI *SZER* (DEPS*  *2 ) ) ) 

A2=CMPLX(0. ,B*DSIG/(4 . *SZEK*FI*DEPS) ) 

A3=1.+2.*A1 
A4=2.*C*A1 
A5=l .+A4 
CJ( 1 )=-2.*C*AJ 
BJ( 1 )=1 ,+CJ ( 1 ) 

B2=1.-2.*A1 
DO  1000  1=2, N 
BJ( I)=B2 

FJ=l.+B/(2.*FLOAT(I-l )) 
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CJ(I)=-A1*FJ 

AJ(l)=Al*(FJ-2.) 

1000  CONTINUE 
C 

C  PERFORM  THE  FINITE- DIFFERENCE  EVALUATION  AS  MANY  TIMES  AS  IS 
C  REQUIRED. 

C 

DO  3500  11=1, NT 
DO  1500  1=1, N 
P(I)=W(I) 

1500  CONTINUE 

P(N+1)=1.5*P(N)-0.5*P(N-1) 

D(1)=W(1)*A5-W(2)*A4 
EPS=0 . 

DO  2000  1=2, N 
EPS=EPS+DEPS 

D( I )=P(I-1 )*(— A1+A2/EPS )+P( I) *A3+P( 1+1 )*(— A1-A2/EPS ) 

2000  CONTINUE 

AF( 1 )=BJ( 1  ) 

SF(1)=D(1) 

DO  2500  1=2  ,N 
Rl=AJ(I)/AF(I-l ) 

AF(I)=BJ (I )-Rl *CJ ( 1-1 ) 

SFU)  =  D(I)+R1*SF(I-1  ) 

2500  CONTINUE 
C 

C  ON  EVALUATING  THE  NEW  TIME  ROW,  RESTORE  THE  SAMPLES  TO  THE  W  ARRAY. 
C 

W(N)=SF(N)/AF(N) 

K=N 

DO  3000  1=2, N 
K=K— 1 

W(K)=(SF(K)+CJ(K)*W(K+1))/AF(K) 

3000  CONTINUE 
3500  CONTINUE 

IF  (RG  .LT.  0.5*RSPH)  GO  TO  8888 
8888  RETURN 
C 

C  SPHERICAL  SPREADING  BY  IMPOSITION  —  FIRST 

C  APPLY  THE  RANGE  DECREMENT  CG,  THEN  INTERPOLATE  TO  ACCOUNT  FOR 
C  THE  RECTANGULAR  SAMPLE  GRID. 

C 

5000  RI=DS/(S2ER*FI) 

CG=CMPLX( 1 . ,RG)/CMPLX( 1 . ,RG+RI) 

DO  5200  1=1, N 
P( I)=W( I)*CG 
5200  CONTINUE 
W(1)=P(1) 

DG=RG/(RG+RI) 

IF  (RG  .LE.  (SS3*RSPH) )  DG=RG/(RG+SS4*RI ) 

IF  (RG  .LE.  (SS1*RSPH))  DG=RG/ ( RG+SS2*RI) 

EG=0 . 

DO  5500  1=2, N 

EC=EG+DC 

J=IFIX(EG)+) 

K=J+1 

FG=EG-FLOAT(J-l) 

HG=1.-FG 

W( I )=P( J)*HG+P(K)*FC 
5500  CONTINUL 

GO  TO  8888 
END 


SUBROUTINE  SSWTCH(N,V) 

LOGICAL  V, S(  13) /.FALSE. ,  .FALSE.,  .FALSE.  ,  .F,\LSE.  ,  .FALSE., 

1  .TRUE.,  .FALSE. , .TRUE. ,  .FALSE., 

2  .FALSE.,. FALSE.,. TRUE.,  .FALSE./ 
L=N+1 

V=. FALSE. 
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IF  <(L  .GT.  13)  .OR.  (L  .LT.  1))  CO  TO  8888 
V=S(L) 

8888  RETURN 
END 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

1000 

c 

c 

c 

c 

c 

3000 


SUBROUTINE  OFT  (  Z,  W,  ISQ,  N,  NC  ) 

DISCRETE  FOURIER  TRANSFORM  SUBROUTINE  WITH  PASSED  ARGUMENTS 

THIS  SUBROUTINE  IS  OPTIMIZED  FOR  RAPID  SEQUENTIAL  PERFORMANCE 
OF  SUCCESSIVE  DIRECT  AND  INVERSE  TRANSFORMS  OF  THE  SAME  SIZE. 


DEFINITIONS  OF  ARGUMENTS: 

Z  -  ON  CALL,  INPUT  COMPLEX  SEQUENCE;  ON  RETURN,  OUTPUT 
COMPLEX  SEQUENCE. 

W  -  COMPLEX  WORKING  ARRAY;  THESE  VALUES  MUST  NOT  BE  CHANGED 
BETWEEN  CALLS  TO  THE  SAME  SIZE  DFT. 

ISQ  -  INTEGER  ARRAY  USED  TO  STORE  THE  SEQUENCE  NUMBERS 
WHICH  ARE  EXCHANGED  IN  THE  BIT  REVERSAL. 

N  -  NUMBER  OF  POINTS  IN  THE  FOURIER  TRANSFORM;  MINIMUM 
LENGTH  OF  THE  Z,  W,  AND  ISQ  ARRAYS. 

NC  -  CONTROL  INTEGER  VARIABLE;  LESS  THAN  ZERO  FOR  INVERSE 
TRANSFORM,  OTHERWISE  DIRECT  IS  ASSUMED. 

THE  MINIMUM  ALLOWABLE  N  IS  8  AND  THE  MAXIMUM  IS  4096;  THE  VALUE 
OF  N  MUST  BE  AN  INTEGER  POWER  OF  2. 


COMPLEX  C,  D,  Z(N),  W(N) 

REAL*8  CO,  SI.  CD,  SD,  SS .  ARG,  SO 
INTEGER  JU(12),  JD(12),  ISQ(N),  OLDN/-9999/ 

DETERMINE  IF  THE  VALUE  OF  N  CALLED  FOR  IS  CURRENT. 

IF  SO,  CHOOSE  DIRECT  TRANSFORM  (IOFS=0)  OR  INVERSE  TRANSFORM. 

IF  (OLDN  .NE.  N)  GO  TO  7000 
IOFS=0 

IF  (NC  .LT.  0)  I0FS=N2 


FIRST  DFT  LOOP;  MULTIPLY  BY  AMPLITUDE  FACTOR  AND  TAKE  THE 
SUM  AND  DIFFERENCE.  SINCE  THE  COMPLEX  EXPONENTIAL  SAMPLE 
FOR  THIS  LOOP  IS  1,  COMPLEX  MULTIPLICATION  IS  NOT  REQUIRED. 


DO  3100  1=1, N2 
II=I+N2 


d=z[ii) 


C=C*AMPL 


D=D*AMPL 


Z(I)=C+D 
Z(II)=C-D 
3100  CONTINUE 


C 

C 

C 

C 


4500 

C 

C 

C 


BIT-REVERSED  RESEQUENCING  OF  THE  Z  ARRAY. 
STORED  IN  ISQ  TO  CONTROL  THE  SEQUENCING. 


IN2=N2+1 
DO  4500  1=2, N2 
IN2=IN2+1 


IF 

C 

Z 


.LT.  0)  GO  TO 


F  (ISQ(I)  . 

=Z( ISO ( I )  ) 
(ISQ(I))=Z(ISQ(IN2)) 
Z( ISQ ( IN2) )=C 
CONTINUE 


4500 


REMAINING  ITERATIONS  OF  THE  DFT. 

DO  5100  1=2, L2N 
K=JU ( I ) 

KK=JD( I ) 


USE  THE  VALUES 


DO  5100  11=1, KK 

M=M+KI 

J=1-KK 

DO  5100  III=1,K 

J=J+KK 

L=III+M 

LL=L+K 

C=Z(L) 

D=Z(LL) 

IF  (J.EQ.l)  GO  TO  5000 
D=D*W( J+IOFS) 

5000  Z(L)=C+D 
Z(LL)=C-D 
5100  CONTINUE 
C 

C  RETURN  AFTER  TRANSFORM  IS  COMPLETED. 

C 

9999  RETURN 
C 

C  EVALUATE  THE  COMPLEX  EXPONENTIAL  SAMPLES  IN  THE  W  ARRAY,  THE 
C  POWERS  OF  2  IN  ASCENDING  AND  DESCENDING  ORDERS  RESPECTIVELY 

C  IN  JU  AND  JD,  AND  THE  BIT-REVERSE  PAIRS  FOR  SWAPPING  IN  ISQ. 

C 

7000  OLDN=N 

SS=DFLOAT(N) 

ARG=- 6 . 283 185  307 1 8D0 

SS=ARG/SS 

SD=DSIN(SS) 

CD=DCOS(SS) 

CO=CD 

SI=SD 

N8=N/8 

N4=N8+N8 

N2=N4+N4 

N4P2=N4+2 

N2P2=N2+2 

W(1)=CMPLX(1. ,0.) 

J=1 

OLDO=0 

W(N4+I)=CMPLX(0. ,-l.) 

7200  J=J+1 

COP=SNGL(CO) 

SIP=SNGL(SI) 

W(J)=CMPLX(COP.SIP) 

W(J+N4)=CMPLX(SIP,-COP) 

IF  (J.GT.N8)  GO  TO  7300 
W(N2P2-J)=CMPLX(-COP, SIP) 

W(N4P2-J)=CMPLX(-SIP,-C0P) 

SO=SI 

SI=SI*CD+CO*SD 
CO=CO*CD-SO*SD 
GO  TO  7200 
C 

C  DETERMINE  THE  TRANSFORM  NORMALIZATION  FACTOR  AMPL  AND  FILL  IN 
C  THE  REMAINDER  OF  THE  W  ARRAY. 

C 

7300  AMPL=l./SQRT(FLOAT(N)) 

K=N2 

DO  7400  1=1, N2 
K=K+1 

W(K)=CONJG(W(I) ) 

7400  CONTINUE 
C 

C  DETERMINE  THE  LOG  TO  THE  BASE  2  OF  N,  L2N ,  AND  FILL  IN  THE  JU 
C  AND  JD  ARRAYS. 

C 

KK=  1 

A=AL0G10( FLOAT( N ) )*3 .321928  +  0.5 
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L2N=IFIX(A) 

LL=L2N 

DO  7500  L=1,L2N 
JU(L)=KK 
JD(LL)=KK 
LL=LL- 1 
KK=KK.+KK 
7500  CONTINUE 
C 

C  FILL  IN  THE  ISQ  ARRAY  WITH  THE  BIT-REVERSE  PAIRS. 
C 

DO  8100  1=1, N 
ISQ(I)=-1 
8100  CONTINUE 
MM=1 

DO  8500  1=2, N 

K=I-1 

KK=0 

DO  8200  L=1,L2N 
JJ=K-JD(L) 

IF  (JJ  .LT.  0)  GO  TO  8200 
KK=KK+JU(L) 

K=JJ 

IF  (K  .LE.  0)  GO  TO  8300 
8200  CONTINUE 
8300  KK=RK.+  1 
JJ=KK-I 

IF  (JJ  .LE.  0)  GO  TO  8500 
MM=MM+1 
ISQ(MM)=I 
ISQ(MM+N2)=KK 
8500  CONTINUE 

GO  TO  1000 
END 


i 

. 

i 
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PROGRAM  PRPGT2 

LARGE  FINITE-AMPLITUDE  WAVE  PROPAGATION  PROGRAM  —  RELEASE  SIX 
AS  REVISED  IN  FEBRUARY  1981. 

DEFINITION  OF  MAJOR  PROGRAM  VARIABLES: 

UZERO  -  CONTAINS  THE  INITIAL  OR  BASE  WAVEFORM,  BEFORE  THE 
NONLINEAR  OPERATOR 

USIGMA  -  CONTAINS  THE  WAVEFORM  AFTER  THE  NONLINEAR  OPERATOR 
ALIST  -  ARRAY  FOR  STORING  LISTS  OF  REAL  NUMBERS  ON  INPUT 
Z  -  COMPLEX  ARRAY  FOR  THE  SPECTRUM 

W  -  COMPLEX  ARRAY  FOR  THE  VALUES  OF  EXP(J2Pl/N) 

ISQ  -  INTEGER  ARRAY  FOR  THE  BIT-REVERSAL  PAIRS. 

CALFA  -  COMPLEX  ARRAY  FOR  THE  FACTORS  DEFINING  THE  ATTENUATION 

AND/OR  DISPERSION  DURING  PROPAGATION 
ATTN  -  .TRUE.  IF  ATTENUATION  IS  ENABLED 

DISP  -  .TRUE.  IF  DISPERSION  IS  ENABLED 

RE  LX  -  .TRUE.  IF  RELAXATION  IS  ENABLED 

AMPLR  -  .TRUE.  IF  THE  OPERATOR  IS  NOT  TO  COMPENSATE  FOR 

FINITE  AMPLITUDE  LOSSES  —  THIS  IS  THE  DEFAULT 
ATRK  -  ARRAY  TO  STORE  SPECTRAL  POINTS  FOR  PRINTED  DISPLAY 

ITRK  -  ARRAY  TO  STORE  FREQUENCIES  OF  THE  SPECTRAL  POINTS 

TO  B5  PRINTED 

ILIST  -  ARRAY  TO  STORE  DFT  CELL  NUMBERS  OF  SPECTRAL  POINTS 
TO  BE  PRINTED 

REVISION  II  JULY  1980  -  AMPLITUDE  NORMALIZATION  TO  A  SPECIFIED  NUMB 
OF  EQUAL  INPUT  COMPONENTS  PROVIDED  VIA  INPUT  VARIABLE 
'IDUM'. 

REAL  UZERO (2048) ,  USIGMA(2048) ,  ALIST(IO),  ATRK (II),  RALFA(2048), 

1  RZ(4096) 

COMPLEX  Z( 2048) ,  W(2048),  CALFA(1024),  CF 
LOGICAL  ATTN.  DISP,  RE LX.  LV,  AMPLR,  S04 ,  S05 
INTEGER  IDAT(2),  ITRK(U)  ,  tLIST(li),  I$Q(2048) 

EQUIVALENCE  (Z( 1 ) ,RZ( 1 )  ),  (CALFA( 1 ) ,RALFA( 1 ) ) 

READY  TO  RUN. 

OO  LPPAGE=0 
FMUL=1. 

CALL  DATE(IDAT) 

READ( 5 , 802 0 , END=7  000 )  ICARD 

INPUT  CASE  NUMBER  ICASE,  NUMBER  OF  FOURIER  SAMPLES  IFOUR,  AND 
SCALED  FOURIER  FUNDAMENTAL  FSUBF . 

000  READ(ICARD,8020,END=7000)  ICASE,  IFOUR,  FSUBF 

IMAXF=IFOUR/2 

INPUT  MEDIUM  PARAMETERS: 

ADELTA  -  ATTENUATION  COEFFICENT; 

CZERO  -  SMALL-SIGNAL  SPEED  OF  SOUND; 

BETA  -  PARAMETER  OF  NONLINEARITY; 

EMM  -  SIMPLE  DISPERSION  PARAMETER; 

TSUBC  -  RELAXATION  TIME  OF  MEDIUM,  TIMES  SCALING  FREQUENCY. 

READ(ICARD,8000,END=7000)  ADELTA,  CZERO,  BETA,  EMM,  TSUBR 

SENSE  SWITCH  OPTIONS: 

12  -  ON  TO  ENABLE  DISPERSION 
11  -  ON  TO  ENABLE  RELAXATION 
10  -  ON  TO  ENABLE  ATTENUATION 

9  -  ON  TO  INPUT  SIGNAL  PARAMETERS,  OFF  FOR  SCALED  PARAMETERS 
EPSLON  -  SIGNAL  MACH  NUMBER 
FREQ  -  SIGNAL  FREQUENCY  IN  HERTZ 
GAMMA  -  SCALED  ACOUSTIC  REYNOLD'S  NUMBER 
CLAMBD  -  SCALED  RELAXATION  PARAMETER 
SMALLD  -  SCALED  SIMPLE  DISPERSION  PARAMETER 
RSUBC  -  PLANE-WAVE  CRITICAL  RANGE  OR  SHOCK 
FORMATION  DT STANCE 


l 
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8  -  INV ISCID  OPERATOR  DISABLED 

7  -  ON  TO  DISABLE  NORMALIZATION  OF  INITIAL  FUNCTION. 

CALL  SSWTCH(12,DISP) 

CALL  SSWTCHi 1 1 , RE LX) 

CALL  SSWTCHC 10, ATTN) 

CALL  SSWTCHC 9, LV) 

IF  (LV)  GO  TO  1500 

ENTER  SCALED  PARAMETERS 

READ(ICARD,8000,END=7000)  SMALLD,  CLAMBD ,  GAMMA 
COMPUTE  THE  CORRESPONDING  SIGNAL  PARAMETERS. 

EPSLON=EMM* CLAMBD/ (2.* BETA) 

IF  (GAMMA  .NE.  0.)  FREQ  =  6 . 28 3I853*BETA*EPSLON  /  (CZERO* 

1  ADELTA*GAMMA) 

IF  (FSUBF  .LT.  0.)  FSUBF  =  FREQ 
RSUBC  =  CZERO/ (6. 283185*FREQ*BETA*EPSLON) 

GO  TO  1600 

INPUT  EPSLON  AND  FREQ. 

READ(ICARD,8020,END=7000)  EPSLON,  FREQ 

COMPUTE  CORRESPONDING  SCALED  PARAMETERS. 

IF  (FSUBF  .LT.  0.)  FSUBF  =  FREQ 
CLAMBD  =2 . *BETA*EPSLON/EMM 

GAMMA=CLAMBD*EMM*3 . 141593/ (CZERO*ADELTA*FSUBF) 

SMALLD  =  CLAMBD/ (6 . 283 1853*TSUBR*FSUBF ) 

RSU BC=CZERO/ ( 6 . 28  31853*FREQ*BETA*EPSLON ) 

SET  UP  FOR  ITERATION 

ISC  -  INDEX  AT  WHICH  SIGMA  CHANGES;  AT  THIS  STEP  NUMBER  THE 

RANGE  STEP  SIZE  IS  MULTIPLIED  BY  SIGST2  (OR  BY  10  IF 
SIGST2  IS  LESS  THAN  0). 

ISF  -  FILE  TO  WHICH  SPECTRA  ARE  TO  BE  WRITTEN;  LESS  THAN  ZERO 
TO  SUPPRESS  SPECTRAL  OUTPUT. 

IPR  -  NUMER  OF  RANGE  STEPS  PER  PRINTOUT  CYCLE. 

IDUM  -  IF  NOT  ZERO,  THE  ASSUMED  NUMBER  OF  EQUAL  AMPLITUDE 

INPUT  COMPONENTS,  USED  FOR  NORMALIZATION  OF  THE  SPECTRAL 
AMPLITUDE  OUTPUT.  IF  ZERO,  NO  NORMALIZATION  IS  PERFORMED. 
SIGSTP  -  INITIAL  NUMBER  OF  RANGE  STEPS  PER  UNIT  SIGMA. 

SIGST2  -  IF  GREATER  THAN  ZERO,  MULTIPLE  OF  RANGE  STEP  TO  BE 
USED  AFTER  STEP  ISC. 

SMAX  -  MAXIMUM  SCALED  RANGE  OF  INTEREST. 

1600  READ(ICARD,8060)  ISC,  ISF,  IPR,  IDUM,  SIGSTP,  SIGST2 ,  SMAX 
AMPNRM=1 . 

IF  (IDUM  .GT.  0)  AMPNRM=FL0AT(IDUM)/SQRT(FL0AT(IF0UR/4)) 

CALL  SSWTCH(8.AMPLR) 

AMPLR= . NOT . AMP  LR 
ATTF=0 . 

ATT2=Q. 

DSPF=0. 

TF2=(TSUBR*FSUBF) **2 
STEP  1=1. /SIGSTP 

IF  (ATTN  .AND.  GAMMA  .NE.  0.)  ATTF=STEP I/GAMMA 
IF  (RELX)  GO  TO  1700 

IF  (DISP  .AND.  SMALLD  .NE.  0.)  DSPF=-STEPI/SMALLD 
GO  TO  1800 

1700  IF  (CLAMBD  .EQ.  0.)  GO  TO  1800 
DSPF  =  STEPI*TF2/CLAMBD 

IF  (ATTN)  ATT 2  =  TSUBR*FSUBF*GAMMA/ CLAMBD 
1800  READ( ICARD.8040)  LI,  ILIST 
NTRK=-J 

IF  (LI  .LT.  1)  CO  TO  2000 
IF  (LI  .GT.  11)  LI=ll 


M 
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"1 


1900 

C 

C 

C 

2000 


C 

C 

C 

2100 


C 

C 

C 

C 

C 

C 

C 


C 

C 

C 

2300 


2400 


C 

c 

c 


2500 


NTRK=LI 

DO  1900  1=1, LI 
ITRK( I )=ILIST( I )+l 
CONTINUE 


CET  INITIAL  WAVEFORM 


LPPAGE=LPPAGE+1 

WRITE(6,9090)  ICASE,  IDAT ,  LPPAGE 
WRITE(b,9010)  ICARD 

CALL  UN IT( IFOUR, UZERO,FSUBF, FREQ, FMUL, GAMMA, CLAMBD , EPSLON , SMALLD 
1  ICARD) 

IF  (UZERO(l)  .LE.  1.)  GO  TO  2100 
WRITE(6 , 9180 ) 

GO  TO  6000 


SET  UP  ITERATION  VARIABLES 

SIGMA  =  0. 

KOUNT=0 

DS IG=1 . /SIGSTP 

DS=FLOAT( IFOUR) *DSIG/6. 283 1853 

IFH  =  IFOUR/ 2 

IF2  =  IFH+1 

IF21  =  IFOUR+2 

KPR=0 

GET  THE  COMPLEX  FACTORS  ACCOUNTING  FOR  ATTENUATION,  DISPERSION, 
AND  RELAXATION. 

CALL  CMP FAC( C ALFA , IFOUR , ATTF , ATT2 , DSPF , TF2 , RELX , IMAXF , ICARD ) 

DISPLAY  RUN  PARAMETERS 

WRITE! 6, 9 100)  DISP,  ATTN,  RELX,  AMPLR 

WRITE! 6,91 10)  GAMMA,  CLAMBD,  SMALLD 

WRITE! 6, 9120 )  EPSLON,  FREQ,  FSUBF 

WRITE! 6,9 130)  DSPF,  ATTF,  ATT2 

WRITE! 6, 9140 )  IFOUR,  DS,  SIGSTP 

WRITE! 6,9150)  BETA,  EMM,  CZERO 

WRITE! 6,9160;  ADELTA,  TSUBR 

IF  (ISC  .GT.  0)  WRITE( 6, 9190)  SIGST2 ,  ISC 

IF  (ISF  .GE.  0)  WRITE! 6,9040)  ISF 

IF  (IPR  .GT.  1)  WRITE! 6, 9050)  IPR 

LPPAGE=LPPAGE+1 

WRITE! 6, 9090)  ICASE,  IDAT,  LPPAGE 

KLINPR=1 

WRITE (6, 9080) 

IF  (NTRK  .GT.  0)  WRITE(6,9200)  (ILIST(KK) ,KK=1 , NTRK) 

COMPUTATION  OF  U(SIGMA)  GIVEN  U(ZERO) 

IF  (KOUNT  .NE.  ISC)  GO  TO  2400 
WRITE! 6, 92 10)  KOUNT,  SIGST2 
DS  =DS  *SIGST2 
ATTF=ATTF*SIGST2 
DSIG=DSIG*SIGST2 
DSPF=DSPF*SIGST2 

CALL  CMP FAC( C ALFA, IFOUR ,ATTF ,ATT2 , DSPF, TF2 , RELX, IMAXF, ICARD) 

KOUNT=KOUNT+l 

SIGMA  =  SIGMA  +  DSIG 

KPR=KPR+1 

IF  (AMPLR)  GO  TO  2600 

GET  RMS  AMPLITUDE  BEFORE  NON  LINEAR  OPERATOR 

RMSBEF  =  0. 

DO  2500  1=) .IFOUR 

RMSBEF  =  RMSBEF  +  UZERO(I)**2 

CONTINUE 


i 


138 


1* 


RMSBEF=SQRT( RMS BEF/FLOAT( IFOUR) ) 

c 

C  FINITE  DIFFERENCE  INTERPOLATOR 
C 

2600  DX=0. 

DO  2800  1=1, IFOUR 
DX=DX+l. 

DL=UZERO(I)*DS 

XT=DX+DL 

M=XT 

XT=XT-FLOAT(M) 

MO=M 

IF(MO  .LT.  1)  MO=MCH-IFOUR 
IF  (MO  .GT.  IFOUR)  MO=MO-IFOUR 
M1=M+1 

IF (XT  .LT.  0.)  M1=M-1 
IF  (MI  .LT.  1)  M1=M1+ IFOUR 
IF  (Ml  .GT.  IFOUR)  M1=M1-IF0UR 
IF  (XT  .LT.  0.)  XT=-XT 

USIGMA(I)  =  UZERO(MO)*(1 .-XT)  +  UZER0(M1)*XT 
2800  CONTINUE 
C 

C  END  OF  THE  F.D.I.  LOOP. 

C 

IF  (AMPLR)  GO  TO  3100 
C 

C  GET  RMS  AMPLITUDE  AFTER  NON  LINEAR  OPERATOR 
C 

RMSAFT  =  0. 

DO  2900  1=1, IFOUR 
RMSAFT=RMSAFT  +  USIGMA(I)**2 
2900  CONTINUE 

RMSAFT  =  SQRT(RMSAFT/FLOAT( IFOUR)) 

C 

C  FIX  UP  THE  DIFFERENCE 

C 

FACTOR  =  RMSBEF/RMSAFT 
DO  3000  1=1, IFOUR 
USIGMA(I)  =  USIGMA(I)  *  FACTOR 
3000  CONTINUE 
C 

C  DISPERSION  AND/OR  ATTENUATION,  IF  EITHER,  ARE  APPLIED  IN  THE 
C  FREQUENCY  DOMAIN. 

C  IF  EITHER  ATTENUATION  OR  DISPERSION  IS  ENABLED,  THEN  THE 

C  SPECTRUM  MODIFICATION  MUST  BE  PERFORMED--  IF  NIETHER  IS  ENABLED, 
C  TRANSFORM  INTO  THE  SPECTRAL  DOMAIN  ONLY  IF  SPECTRAL  INFORMATION 
C  IS  TO  BE  OUTPUT. 

C 

3100  IF  (ATTN  .OR.  DISP)  GO  TO  3200 
IF  (NTRK  .LT.  1 )  GO  TO  5200 
IF  (RPR  .NE.  IPR)  GO  TO  5200 
C 

C  UNPACK  THE  REAL  DATA  INTO  THE  REAL  PARTS  OF  A  COMPLEX  ARRAY. 

C 

3200  DO  3300  1=1, IFOUR 

Z(I)  =  CMPLX(USIGMA( I) ,0.) 

3300  CONTINUE 

CALL  DFT  (  Z,  W,  ISQ,  IFOUR,  2  ) 

C 

C  BYPASS  THIS  SECTION  IF  ONLY  THE  SPECTRAL  RECORD  IS  DESIRED —  NO 
C  ATTENUATION  OR  DISPERSION. 

C 

IF  (.NOT.  (ATTN  .OR.  DISP)  )  GO  TO  4000 
K=IFOUR 

IF  (DISP)  GO  TO  3500 
C 

C  APPLY  ATTENUATION  ONLY. 

C 

N=1 

DO  3400  1=2 , IMAXF 
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N=N+2 

RN  =  RALFA(N) 

RZ(N)  =  RZ(N)  *  RN 
RZ(N+1)  =  RZCN+l)  *  RN 
Z(K)  =  CONJG(Z(I)) 

K=K~  1 

3400  CONTINUE 

GO  TO  3700 
C 

C  APPLY  DISPERSION  AND  POSSIBLY  ATTENUATION. 

C 

3500  DO  3600  I=2,IMAXF 

Z(I)  =  Z(I)*CALFA(I) 

Z(K)  =  CONJG(Z(I) ) 

K=K— I 

3600  CONTINUE 
C 

C  DELETE  H-F  PART  OF  SPECTRUM 

C 

3700  CF=CMPLX(0. ,0.) 

N=IMAXF+1 
M=IFOUR+2-N 
DO  3800  I=N,M 
Z(I)=CF 
3800  CONTINUE 
C 

C  OUTPUT  SPECTRUM  IF  DESIRED 

C 

4000  IF  (ISF  .LT.  0)  GO  TO  4100 
MMU=  -I 

WRITE  (ISF)  KOUNT,  SIGMA,  IFOUR,  ICASE,  ( I DAT ( KKK ) , KKK= 1,2),  MM  U 
WRITE(ISF)  (Z(KK),KK=1,IF2) 

C 

C  PRINT  TRACKED  SPECTRUM  IF  DESIRED 

C 

4100  IF  (NTRK  . LT.  0)  GO  TO  4300 
DO  4200  1=1, NTRK 

ATRK(I)  =  CABS(Z(ITRK(I)))*AMPNRM 
4200  CONTINUE 
C 

C  INVERSE  FOURIER  TRANSFORM  AND  RESTORE  REAL  FUNCTION  TO  THE 
C  UZERO  ARRAY. 

C 

4300  CALL  DFT  (  Z,  W,  ISQ,  IFOUR,  -2  ) 

J=IFOUR 

DO  4400  1=1, IFOUR 
CF=Z( J) 

F=CABS(CF) 

IF  (REAL(CF)  .LT.  0.)  F=-F 

UZERO(J)  =  F 

J=J-1 

4400  CONTINUE 
C 

C  SEE  IF  A  PRINTED  DATA  RECORD  IS  DESIRED  —  IF  SO,  OUTPUT  IT  AFTER 

C  CHECKING  WHETHER  THE  PAGE  IS  FULL. 

C 

4500  IF  (KPR  .LT.  IPR)  GO  TO  4800 
KPR=0 

KLINP R=KLINPR+1 

IF  (KLINPR  .LE.  50)  GO  TO  4600 

KLINPR=1 

LPPAGE=LPPACE+J 

WRITE(6,9U90)  ICASE,  IDAT,  LPPAGE 
WRITE(6,9080) 

IF  (NTRK  .GT.  0)  WRITE(6,9200)  ( ILIST(KK) ,KK=1 , NTRK) 

WRITE(6,9i 70) 

C 

C  COMPUTE  RMS,  FIND  MAX  AND  MIN,  AND  ANYTHING  ELSE  TO  BE  PRINTED. 

C 

4600  URMS=UZERO( 1)**2 
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4700 


C 

C 

C 

4800 

C 

C 

C 

5000 


C 

C 

C 

5200 

5300 

C 

C 

C 

6000 


C 

C 

C 

7000 


C 

C 

C 

8000 

8020 

8040 

8060 

9000 

9010 

9040 

9050 

9080 

9090 

9100 

9110 

9120 

9130 

9140 


UMIN=UZER0(1) 

UMAX=UMIN 

DO  4700  1=2 , IFOUR 

UI=UZERO(I) 

URMS=  URMS  +  UI**2 
IF  (UI  .GT.  UMAX)  UMAX=UI 
IF  (UI  .LT.  UMIN)  UMIN=UI 
CONTINUE 

URMS=SQRT(URMS/FLOAT( IFOUR)) 

IF  (NTRK  .LT.  1)  WRITE(6,9000)  KOUNT,  SIGMA,  URMS,  UMAX,  UMIN 
IF  (NTRK  .GT.  0)  WRITE(6,9000)  KOUNT,  SIGMA,  URMS,  UMAX,  UMIN, 

1  (ATRK(KK) ,KK=1,NTRK) 

SEE  IF  DESIRED  RANGE  IS  ATTAINED. 

IF  (SIGMA  .GE.  SMAX)  GO  TO  5000 
GO  TO  2300 

END-OF-SEQUENCE  HANDLER 

WRITE(6,9340)  KOUNT,  SIGMA 
KOUNT=- I 

IF  (ISF  .GE.  0)  WRITE  (ISF)  KOUNT,  SIGMA,  IFOUR,  ICASE,  (IDAT(KKK) 
1 , KKK= 1,2).  MMU 
GO  TO  1000 

PUT  NEW  WAVEFORM  BACK  IN  OLD  ARRAY 

DO  5300  1=1, IFOUR 
UZERO(I)  =  US IGMA(I) 

CONTINUE 
GO  TO  4500 

END-OF-STEP  HANDLER 

KOUNT=-l 

IF  (ISF  .GE.  0)  WRITE  (ISF)  KOUNT,  SIGMA,  IFOUR,  ICASE,  (IDAT(KKK) 
l.KKK=l,2),  MMU 
STOP 

ERROR  HANDLER 


KOUNT=-l 

IF  (ISF  .GE.  0)  WRITE  (ISF)  KOUNT,  SIGMA,  IFOUR,  ICASE,  (IDAT(KKK) 
1 , KKK=1 , 2 ) ,  MMU 
STOP 


FORMAT  STATEMENTS 


FORMAT! 8F 10. 1 ) 

FORMAT! 215. 7 F10. 1 ) 

FORMAT! 1615) 

FORMAT14I5 ,6F10. 1) 

FORMAT!'  I5,F8.3,F10.4,2F10.3,11F8.4) 

FORMAT ('OCARD- IMAGE  CONTROL  INPUT  FROM  DEVICE  ,13) 

FORMAT! 'ORECORD  SPECTRA  ON', 13) 

FORMAT! 'ODISPLAY  PRINTED  RECORD  EVERY', 13,'  RANGE  STEPS.') 

FORMAT! '0  STEP  SIGMA  U(RMS)  (MAX)  (MIN)') 

FORMAT! '1  BELLMANS  METHOD  -  RELEASE  FIVE,  7  AUGU 
1ST  1980  -  CASE', 110, 4X,2A4,5X, 'PAGE', 14) 

FORMAT('0  DISPERSION  'Ll.  ;  ATTENUATION  ' ,L1,';  RELAXATION  ' , LI , 

'  INVISCID  OPERATOR  , LI ) 


1 


FORMAT! '0 GAMMA 
1E13.3) 

FORMAT( 'OEPSLON 
1E13.3) 

FORMAT( 'ODSPF 
1E13.3) 

FORMAT( ' OIFOUR 

1) 


vni  >  ua  / 

' , 1PE13.3,' 

' , 1PE13.3,' 

' , 1PE13 . 3, ' 
',113,'  DS 


CLAMBDA  ' ,E13.3,' 
FREQ  '  ,E13 . 3, ' 

ATTF  ' ,E13.3,' 

' , 1PE13.3,' 


SMALL  D  , 

F  SUB  F  ', 

ATT 2  '  , 

SIGSTP  ' , E 13 . 3 


mm mm 


CZERO 


9150 

9160 

9170 

9180 

9190 

9200 

9210 

9340 


FO  RMAT ( ' 0  BETA  1PE13.3,'  EMM  '.E13.3,' 

1E13.3) 

FORMAT! 'ODELTA  '.1PE13.3,'  T  SUB  R  '.E13.3) 

FORMAT! '  ') 

FORMAT! '0***  BEWARE  ***  UNSPECIFIED  INITIAL  WAVEFORM.') 
FORMAT! 'ONOTE:  DELTA  SIGMA  WILL  BE  MULTIPLIED  BY',1PE11.3, 
1  '  AFTER' .16.'  STEPS.') 

FORMAT! 4&X, 11 18) 

FORMAT! 'OAFTER  STEP  ,14,'  MULTIPLY  DELTA  SIGMA  BY',F8.1) 
FORMAT! 'OEND  OF  TASK  AT  STEP', 15,'  SIGMA' , 1PE12 . 3) 

END 


SUBROUTINE  CMPFAC  !  CALFA,  IFOUR,  ATTF,  ATT2 ,  DSPF,  TF2 ,  RE LX, 

1  IMAXF,  ICARD  ) 

C 

C  ATTENUATION,  RELAXATION,  AND  DISPERSION  FACTORS.  SEE  19  X  '77  AND 
C  24  VI  '80  NOTES  FOR  DETAILS. 

C 

C  REVISION  8  VII  1980  --  PROGRAM  MODIFIED  TO  LIST  STYLE  OF 

C  ATTENUATION,  RELAXATION,  AND  DISPERSION  ON  PRINTED  LISTING. 

C 

C  REVISION  30  JULY  1980  —  INPUT  TAKEN  FROM  CARD- IMAGE  FILE  ' IFILE ' . 

C 

COMPLEX  CALFA! IFOUR),  CF 
LOGICAL  RELX 
C 

WRITE! 6,9000) 

CALFA! 1)=1. 

ENN=0 . 

CF=CMPLX(0.  ,0.) 

Nl=IF0UR/2  +  1 
DO  1000  1=2, N1 
ENN=ENN+1. 

EN2=ENN**2 

DEN=1.+EN2*TF2 

EF=EXP(~EN2*ATTF*U  .+ATT2/DEN)  ) 

IF  !EF  .LT.  I.E-IO)  GO  TO  800 

PF=ENN**3*DSPF 

IF  (RELX)  PF=PF/DEN 

CALFA( I)=CMPLX(EF*COS(PF) ,EF*SIN(PF) ) 

GO  TO  1000 
800  CALFA(I)=CF 
1000  CONTINUE 
8888  RETURN 

90O0  FORMAT! 'OCOMPLEX  ATTENUATION,  DISPERSION,  AND  ATTENUATION;  SEE  NOT 
1ES  OF  24  JUNE  1980  FOR  DETAILS.') 

END 


SUBROUTINE  UNIT  (  IFOUR,  UZERO,  FSUBF,  FREQ,  FMUL ,  GAMMA, 

1  CLAMBD ,  EPSLON,  SMALL D  ,  ICARD  ) 

C 

C  WAVEFORM  INITIALIZATION  SUBROUTINE. 

C  DESIGNED/CODED  24  VI  1980. 

C 

C  SOLI  DEO  GLORIA  24  VI  1980  —  F  S  MCKENDREE. 

C 

C  REVISION  30  JULY  1980  —  INPUT  TAKEN  FROM  CARD-IMAGE  FILE  'IFILE'. 
C 

C  REVISION  5  AUGUST  1980  -  ABILITY  TO  DISABLE  INITIAL  FUNCTION 
C  NORMALIZATION  ADUED;  USE  SWITCIi  7  ON  TO  DISABLE. 

C 

REAL  UZERO! IFOUR) 

LOGICAL  SO/ 

C 

C  SET  UP  AND  ZERO  DATA  ARRAY. 

C 

DO  500  1=1, IFOUR 
UZERO! I )=U . 

500  CONTINUE 
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C 

C 

C 

C 

1000 


1500 

C 

C 

C 

2000 


2500 


3000 

8888 

8000 

9000 

9020 

9040 


DP=6 . 2831853*FMUL/FLOAT( IFOUR) 

WRITE(6,9000)  IFOUR,  FSUBF,  FMUL 

READ  FREQUENCY,  AMPLITUDE,  AND  PHASE  OF  EACH  COMPONENT,  AND 
ACCUMULATE  IN  UZERO  ARRAY. 

PI,  F2 ,  A 2,  P2 


P=-DP 

DO  1500  1=1, IFOUR 
P=P+DP 

UZERO ( I ) =UZE  RO ( I )+Al *C0  S ( F 1 *P+Rl )+A2  *COS ( F2  *P+R2 ) 

CONTINUE 

IF  (FI  .GT.  0.  .AND.  F2  .GT.  0.)  GO  TO  1000 

NORMALIZE  TO  UNITY  PEAK  AMPLITUDE. 

CALL  SSWTCH(7,S07 ) 

IF  (S07)  WRITE(6,9040) 

IF  (S07)  GO  TO  8888 
AMAX=ABS(UZERO(l)) 

DO  2500  1=2, IFOUR 
AM=ABS(UZERO(I) ) 

IF  (AM  .GT.  AMAX)  AMAX=AM 

CONTINUE 

AM=1 ./AMAX 

DO  3000  1=1, IFOUR 

UZERO ( I )=UZERO (I ) *AM 

CONTINUE 

RETURN 

FORMAT(8F 10.1) 

FORMAT! 'OUNIT  SUBROUTINE:  IFOUR' ,15,'  FSUBF,  FMUL' . 1P2E12 .3) 
FORMAT! ' 0F1 ,  Al ,  Pi : '  1P3E12 . 3, ' ;  F2,  A2 ,  P2:',3El2.3) 

FORMAT! '0 INITIAL  FUNCTION  WILL  NOT  BE  NORMALIZED  TO  UNITY  PEAK  AMP 
1LITUDE. ' ) 

END 


READ(ICARD ,8000,END=2000)  FI,  Al , 
WRITE( 6 , 9020)  FI,  Al ,  Pi,  F2,  A2, 
Rl=Pl*0l01745329 
R2=P2*0. 01745329 


SUBROUTINE  DFT  (  Z,  W,  ISQ,  N,  NC  ) 

C 

C  DISCRETE  FOURIER  TRANSFORM  SUBROUTINE  WITH  PASSED  ARGUMENTS 

C 

C  THIS  SUBROUTINE  IS  OPTIMIZED  FOR  RAPID  SEQUENTIAL  PERFORMANCE 
C  OF  SUCCESSIVE  DIRECT  AND  INVERSE  TRANSFORMS  OF  THE  SAME  SIZE. 

C 

C  DEFINITIONS  OF  ARGUMENTS: 

C  Z  -  ON  CALL,  INPUT  COMPLEX  SEQUENCE;  ON  RETURN,  OUTPUT 

C  COMPLEX  SEQUENCE. 

C  W  -  COMPLEX  WORKING  ARRAY;  THESE  VALUES  MUST  NOT  BE  CHANGED 

C  BETWEEN  CALLS  TO  THE  SAME  SIZE  DFT. 

C  ISQ  -  INTEGER  ARRAY  USED  TO  STORE  THE  SEQUENCE  NUMBERS 

C  WHICH  ARE  EXCHANGED  IN  THE  BIT  REVERSAL. 

C  N  -  NUMBER  OF  POINTS  IN  THE  FOURIER  TRANSFORM;  MINIMUM 

C  LENGTH  OF  THE  Z,  W,  AND  ISQ  ARRAYS. 

C  NC  -  CONTROL  INTEGER  VARIABLE;  LESS  THAN  ZERO  FOR  INVERSE 

C  TRANSFORM,  OTHERWISE  DIRECT  IS  ASSUMED. 

C 

C  THE  MINIMUM  ALLOWABLE  N  IS  8  AND  THE  MAXIMUM  IS  4096;  THE  VALUE 

C  OF  N  MUST  BE  AN  INTEGER  POWER  OF  2. 

C 

COMPLEX  C,  D,  Z(N),  W(N) 

REAL*8  CO,  SI,  CD,  SD,  SS ,  ARG,  SO 
INTEGER  JU(12),  JD(12),  ISQ(N),  OLDN/-9999/ 

C 

C  DETERMINE  IF  THE  VALUE  OF  N  CALLED  FOR  IS  CURRENT. 

C  IF  SO,  CHOOSE  DIRECT  TRANSFORM  (IOFS=U)  OR  INVERSE  TRANSFORM. 

C 

IF  (OLDN  .NE.  N)  GO  TO  7000 
1000  I0FS=0 

IF  (NC  .LT.  0)  IOFS=N2 
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C 

C 

C 

C 

C 

3000 


FIRST  DFT  LOOP;  MULTIPLY  BY  AMPLITUDE  FACTOR  AND  TAKE  THE 
SUM  AND  DIFFERENCE.  SINCE  THE  COMPLEX  EXPONENTIAL  SAMPLE 
FOR  THIS  LOOP  IS  1,  COMPLEX  MULTIPLICATION  IS  NOT  REQUIRED. 


DO  3100  1=1, N2 
II=I+N2 

c=z(i) 

D=Z(II) 

C=C*AMPL 

D=D*AMPL 


Z(I)=C+D 
Z(II)=C-D 
3100  CONTINUE 


C 

C 

C 

C 


4500 

C 

C 

C 


BIT-REVERSED  RESEQUENCING  OF  THE  Z  ARRAY. 
STORED  IN  ISQ  TO  CONTROL  THE  SEQUENCING. 


IN2=N2+1 

DO  4500  1=2, N2 

IN2=IN2+1 

IF  (ISQ(I)  .LT.  0)  GO  TO 
C=Z(IS^(I)) 


Z(ISQ(I) )=Z(ISQ(IN2) ) 

Z(ISQ(IN2))=C 

CONTINUE 


4500 


REMAINING  ITERATIONS  OF  THE  DFT. 


DO  5100  1=2 
K=JU(I) 

KK= JD ( I ) 
KI=K+K 


,L2N 


USE  THE  VALUES 


M=-KI 

DO  5100  11=1, KK 

M=M+KI 

J=1-KK 

DO  5100  111=1 , K 

J=J+KK 

L=III+M 


LL=L+K 

C=Z(L) 

D=Z(LL) 

IF  (J.EQ.l)  GO  TO  5000 
D=D*W( J+IOFS) 

5000  Z(L)=C+D 
Z(LL)=C-D 
5100  CONTINUE 


C 

C  RETURN  AFTER  TRANSFORM  IS  COMPLETED. 

C 

9999  RETURN 

C 

C  EVALUATE  THE  COMPLEX  EXPONENTIAL  SAMPLES  IN  THE  W  ARRAY,  THE 

C  POWERS  OF  2  IN  ASCENDING  AND  DESCENDING  ORDERS  RESPECTIVELY 

C  IN  JU  AND  JD,  AND  THE  BIT-REVERSE  PAIRS  FOR  SWAPPING  IN  ISQ. 

C 

7000  OLDN=N 

S  S=  DF  LOAT ( N ) 

ARG=-6. 283 18530718D0 

SS=ARG/SS 

SD=DSIN(SS) 

CD=DCOS(SS) 

CO=CD 

SI=SD 

N8=N/8 

N4=N8+N8 

N2=N4+N4 

N4P2=N4+2 

N2P2=N2+2 


1 
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W( 1 )=CMPLX( 1 . ,0 • ) 

J=1 

OLDO=0 

W(N4+1)=CMPLX(0. ,-l.) 

7200  J-J+l 

C0P=SNGL(C0 ) 

SIP=SNGL(S1) 

W( J )=CMPLX(COP ,S IP) 

W(J+N4)=CMPLX(SIP,-COP) 

IF  (J.GT.N8)  GO  TO  7300 
W(N2P2-J)=CMPLX(-COP,SIP) 

W(N4P2-J)=CMPLX(-SIP ,-COP ) 

SO=SI 

SI=SI*CD+CO*SD 
CO=CO*CD-SO*SD 
GO  TO  7200 
C 

C  DETERMINE  THE  TRANSFORM  NORMALIZATION  FACTOR  AMPL  AND  FILL  IN 
C  THE  REMAINDER  OF  THE  W  ARRAY. 

C 

7300  AMPL=l./SQRT(FLOAT(N)) 

K=N2 

DO  7400  1=1, N2 
K=K+1 

W( K)=CONJ  G(W( I ) ) 

7400  CONTINUE 
C 

C  DETERMINE  THE  LOG  TO  THE  BASE  2  OF  N,  L2N,  AND  FILL  IN  THE  JU 
C  AND  JD  ARRAYS. 

C 

KK=1 

A=AL0G10(FL0AT(N) )*3. 321928  +  0.5 
L2N=IFIX(A) 

LL=L2N 

DO  7500  L=1 ,L2N 
JU(L)=KK 
JD(LL)=KK 
LL=LL-1 
KK=KK+KK 
7500  CONTINUE 
C 

C  FILL  IN  THE  ISQ  ARRAY  WITH  THE  BIT-REVERSE  PAIRS. 

C 

DO  8100  1=1, N 
ISQ(I)=-1 
8100  CONTINUE 
MM=1 

DO  8500  1=2, N 

K=I-1 

KK=0 

DO  8200  L=1,L2N 
JJ=K-JD(L) 

IF  (JJ  .LT.  0)  GO  TO  8200 
KK=KK+JU(L) 

K=JJ 

IF  (X  .LE.  0)  GO  TO  8300 
8200  CONTINUE 
8300  KK=KK+1 
JJ=KK-I 

IF  (JJ  .LE.  0)  GO  TO  8500 
MM=MM+1 
ISQ(MM)=I 
IS0(MM+N2)=XK 
8500  CONTINUE 

GO  TO  1000 
END 


SUBROUTINE  SSWTCH(N,V) 

LOGICAL  V,S( 13)/. FALSE., .FALSE. , .FALSE. , .FALSE. , .FALSE. , 

1  .FALSE., .FALSE. , .FALSE., .FALSE. , 

2  .FALSE., .TRUE.  , .FALSE. , .FALSE. / 
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